{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as  plt\n",
    "import pandas as pd\n",
    "from math import log,sqrt,exp\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#根据公式计算期权价值\n",
    "def bsm_call_value(s0,k,t,r,sigma):\n",
    "    d1 = ( log( s0/k ) + ( r + 0.5*sigma**2 )*t )/( sigma*sqrt(t) )\n",
    "    d2 = ( log( s0/k ) + ( r - 0.5*sigma**2 )*t )/( sigma*sqrt(t) )\n",
    "    value = ( s0*stats.norm.cdf( d1,0.,1. ) - k*exp( -r*t )*stats.norm.cdf( d2,0.,1 ))\n",
    "    #print('cvalue',value)\n",
    "    return value\n",
    "\n",
    "#求vega\n",
    "def bsm_vega(s0,k,t,r,sigma):\n",
    "    d1 = log( s0/k ) + ( r + 0.5*sigma**2 )*t/( sigma*sqrt(t) )\n",
    "    vega = s0*stats.norm.cdf(d1,0.,1.)*sqrt(t)\n",
    "    #print('vega',vega)\n",
    "    return vega\n",
    "\n",
    "#牛顿迭代法求隐含波动率，迭代次数设为100\n",
    "def bsm_call_imp_vol_newton(s0, k, t, r, c0, sigma_est, it = 100):    \n",
    "    for i in range(it):\n",
    "        sigma_est -= ((bsm_call_value(s0, k, t, r, sigma_est) - c0)/ \n",
    "                      bsm_vega(s0, k, t, r, sigma_est))\n",
    "    return sigma_est"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#二分法求隐含波动率\n",
    "def bsm_call_imp_vol_dichotomy(s0,k,t,r,c):\n",
    "    c_est = 0\n",
    "    top = 3  #波动率上限\n",
    "    floor = 0  #波动率下限\n",
    "    sigma = ( floor + top )/2 #波动率初始值\n",
    "    \n",
    "    while abs( c - c_est ) > 1e-8:\n",
    "        c_est = bsm_call_value(s0,k,t,r,sigma) \n",
    "        #根据价格判断波动率是被低估还是高估，并对波动率做修正\n",
    "        if c - c_est > 0: #f(x)>0\n",
    "            floor = sigma\n",
    "            sigma = ( sigma + top )/2\n",
    "        else:\n",
    "            top = sigma\n",
    "            sigma = ( sigma + floor )/2\n",
    "    return sigma \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "#读取行情数据\n",
    "call_50etf_16 = pd.read_csv('D:\\\\data\\\\options\\\\call_50etf_16.csv',encoding='utf-8')\n",
    "call_50etf_51 = pd.read_csv('D:\\\\data\\\\options\\\\call_50etf_51.csv',encoding='utf-8')\n",
    "call_50etf_79 = pd.read_csv('D:\\\\data\\\\options\\\\call_50etf_79.csv',encoding='utf-8')\n",
    "\n",
    "#取期权价格和行权价\n",
    "price_16 = call_50etf_16['最新价']\n",
    "k_16 = call_50etf_16['行权价']\n",
    "price_51 = call_50etf_51['最新价']\n",
    "k_51 = call_50etf_51['行权价']\n",
    "price_79 = call_50etf_79['最新价']\n",
    "k_79 = call_50etf_79['行权价']\n",
    "\n",
    "#到期时间年化\n",
    "t_16 = 16/365\n",
    "t_51 = 51/365\n",
    "t_79 = 79/365\n",
    "\n",
    "#标的初始价格\n",
    "s0 = 2.533\n",
    "#无风险利率，用shibor近似\n",
    "rf = 0.025\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "imp_vol_newton_16:\n",
      "[0.3638689603157307, 0.3647065048442553, 0.33796965226626546, 0.3087508280418194, 0.29010352112059345, 0.27704783086782625, 0.2678504434116786, 0.27958298121452063, 0.2815159982582914, 0.2892477867510613, 0.2948778490238011, 0.3001603015476018, 0.3079952199130588]\n",
      "imp_vol_dichotomy_16:\n",
      "[0.3638690114021301, 0.3647065758705139, 0.33796969056129456, 0.3087505102157593, 0.290103480219841, 0.2770475149154663, 0.2678511142730713, 0.2795829623937607, 0.28151603043079376, 0.2892477214336395, 0.29487812519073486, 0.3001604676246643, 0.30799537897109985]\n",
      "-------------------------------------------------------------------------------------------------------------------------------------\n",
      "imp_vol_newton_51:\n",
      "[0.2663830055838, 0.2623349770104388, 0.25886692974504383, 0.25514765329363565, 0.25388257585478174, 0.25112748137747243, 0.25311127192458244, 0.25282431711916714, 0.25334527303228827]\n",
      "imp_vol_dichotomy_51:\n",
      "[0.26638298481702805, 0.26233501732349396, 0.25886694341897964, 0.2551477253437042, 0.25388259440660477, 0.2511274963617325, 0.25311121344566345, 0.2528250217437744, 0.25334523618221283]\n",
      "-------------------------------------------------------------------------------------------------------------------------------------\n",
      "imp_vol_newton_79:\n",
      "[0.30525607059637466, 0.2900199629942286, 0.29426758619557364, 0.26917014044939297, 0.2708636644397624, 0.25993171220348904, 0.2610946078237793, 0.2527101923035752, 0.2523987108427234, 0.24752459839396573, 0.2483057133362068, 0.24353321749378323, 0.24308705976680656, 0.23851008170784155, 0.2406859629151288, 0.24016117575185575]\n",
      "imp_vol_dichotomy_79:\n",
      "[0.30525608360767365, 0.29002001881599426, 0.29426760971546173, 0.2691701799631119, 0.2708636373281479, 0.25993166863918304, 0.261094618588686, 0.2527102008461952, 0.25239837169647217, 0.2475244402885437, 0.24830570071935654, 0.2435331791639328, 0.243087038397789, 0.2385101616382599, 0.24068592488765717, 0.24016119539737701]\n",
      "-------------------------------------------------------------------------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "#用两种方法分别求三只期权的隐含波动率，并打印\n",
    "sigma_init=1\n",
    "sigma_16_newton=[]\n",
    "sigma_16_dichotomy=[]\n",
    "for i in range(call_50etf_16.shape[0]):\n",
    "    sigma_16_newton.append(bsm_call_imp_vol_newton(s0,k_16[i],t_16,rf,price_16[i],sigma_init))\n",
    "    sigma_16_dichotomy.append(bsm_call_imp_vol_dichotomy(s0,k_16[i],t_16,rf,price_16[i]))\n",
    "    \n",
    "print('imp_vol_newton_16:')\n",
    "print(sigma_16_newton)\n",
    "print('imp_vol_dichotomy_16:')\n",
    "print(sigma_16_dichotomy)\n",
    "print('-------------------------------------------------------------------------------------------------------------------------------------')\n",
    "#\n",
    "sigma_51_newton=[]\n",
    "sigma_51_dichotomy=[]\n",
    "for i in range(call_50etf_51.shape[0]):\n",
    "    sigma_51_newton.append(bsm_call_imp_vol_newton(s0,k_51[i],t_51,rf,price_51[i],sigma_init))\n",
    "    sigma_51_dichotomy.append(bsm_call_imp_vol_dichotomy(s0,k_51[i],t_51,rf,price_51[i]))\n",
    "\n",
    "print('imp_vol_newton_51:')    \n",
    "print(sigma_51_newton)\n",
    "print('imp_vol_dichotomy_51:')\n",
    "print(sigma_51_dichotomy)\n",
    "print('-------------------------------------------------------------------------------------------------------------------------------------')\n",
    "\n",
    "#\n",
    "sigma_79_newton=[]\n",
    "sigma_79_dichotomy=[]\n",
    "for i in range(call_50etf_79.shape[0]):\n",
    "    sigma_79_newton.append(bsm_call_imp_vol_newton(s0,k_79[i],t_79,rf,price_79[i],sigma_init))\n",
    "    sigma_79_dichotomy.append(bsm_call_imp_vol_dichotomy(s0,k_79[i],t_79,rf,price_79[i]))\n",
    "\n",
    "print('imp_vol_newton_79:')     \n",
    "print(sigma_79_newton)\n",
    "print('imp_vol_dichotomy_79:')\n",
    "print(sigma_79_dichotomy)\n",
    "print('-------------------------------------------------------------------------------------------------------------------------------------')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEKCAYAAADjDHn2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XlcVNX/x/HXmYFhB3dcUERw38Vt3ADJLTPbtFxLs1XbbPXXNyvbM9NKy3I3UbOstNTMEFwCTVHcRVFR0cxERRBhtvP7Y0ZFQ0BgmAHO8/G4D2bu3HvnPZTz4dx7zzlCSomiKIqi5Efj6ACKoiiK81PFQlEURSmQKhaKoihKgVSxUBRFUQqkioWiKIpSIFUsFEVRlAKpYqEoiqIUSBULRVEUpUCqWCiKoigFcnF0gJJSrVo1Wb9+/SLvf/nyZby8vEouUAlz9nzg/BmdPR+ojCXB2fOBc2VMSEg4J6WsXuCGUspysYSGhsriiImJKdb+9ubs+aR0/ozOnk9KlbEkOHs+KZ0rI7BdFuI7Vp2GUhRFUQqkioWiKIpSIFUsFEVRlAKVmwvciqJUPEajkdTUVLKzs6+t8/Pz48CBAw5MVTBHZHR3dycgIABXV9ci7a+KhaIoZVZqaio+Pj7Ur18fIQQAGRkZ+Pj4ODhZ/ko7o5SStLQ0UlNTCQoKKtIx1GkoRVHKrOzsbKpWrXqtUCh5E0JQtWrVG1pgt0sVi7IgPp56UVEQH+/oJIridFShKJzi/p7UaSgnYbZILl0xcinbSPoV63LpignXv7YQMW4o9Y0GiIqC6GjQ6x0dV1GUCkYVCzs4vWY9lpgY/m7bmZONW1/74k/PVQwu2QpCRrZ1fWaOKc9jPR3/Mz0NBjTSgjQYELGxqlgoihMZPXo0v/76KzVq1GDv3r3X1n/xxRdMnz4dFxcX+vfvz8cff5zvcVJSUrjrrrtuOIYzUcWipMXHU3XgnWhNRqpqXfjgoffYUacpAN5uLvi6u+Dr4Yqvhyt1q3ji6+6Kn4crvh4u+HnYHru74udp/Vkt0gvtwO8xZedg0bqgCw937OdTFOUGjzzyCOPGjWPkyJHX1sXExLBixQp2796Nm5sbZ8+edWDCkqGKRUmLjUVnNiGkBa00MyfwMnJCL3zdXXDRFuESUc1wiI7mt8mzmOvVjHcDm9OsxEMrilJUPXr0ICUl5YZ1X331Fa+99hpubm4A1KhRI899ExISGD16NJ6ennTr1u3a+pSUFEaMGMHly5cBmD59Ol26dGHEiBE88MADDBw4EIBhw4bx4IMPEhwczKhRozAYDFgsFpYvX07Dhg1L9HOqYlHSwsMRbjowGBA6HZX79wYvXfGOqdcjnsrm8GYDM2KTmTG0XclkVZRy5O1f9rH/9CXMZjNarbZEjtmsti9vDmh+2/sdOnSITZs28frrr+Pu7s4nn3xChw4d/rPdqFGj+OKLLwgLC+Pll1++tr5GjRqsW7cOd3d3Dh8+zJAhQ9i+fTtjxoxh6tSpDBw4kPT0dOLi4liwYAEvvPACzz33HMOGDcNgMGA2m4v1ufNi17uhhBB9hRBJQohkIcRrebz+pBBijxAiUQixWQjRLNdrrYQQ8UKIfbZt3O2ZtcTo9daL0O+8U6IXo71cBSO7BLJ6z98kn80skWMqimIfJpOJCxcusGXLFiZPnszgwYOxjtl3XXp6OhcvXiQsLAyAESNGXHvNaDTy2GOP0bJlSwYNGsT+/fsBCAsLIzk5mbNnz7JkyRLuv/9+XFxc0Ov1vP/++3z00UccP34cDw+PEv9MdmtZCCG0wAygF5AKbBNCrJRS7s+12WIp5Uzb9ncDnwJ9hRAuwCJghJRylxCiKmC0V9YSp9fb5SL06K5BzN2cwpexyXw6uE2JH19RyrKrLQBn6JQXEBDAfffdhxCCjh07otFoOHfuHK+88go7d+6kRo0aLFu27Ja3s06dOhV/f3927dqFxWLB3f3638ojRowgKiqKpUuXMnfuXACGDh1Kp06dWLVqFX369GH27Nn07NmzRD+TPVsWHYFkKeVRKaUBWAoMzL2BlPJSrqdewNXS2xvYLaXcZdsuTUpZ8u2qMqaqtxvDOtVjReJpTqRlOTqOoii3cM8997B+/XrAekrKYDBQrVo15s2bR2JiIsuXL6dSpUr4+fmxefNmAKKioq7tn56eTq1atdBoNHz77bc3nFZ65JFHmDZtGgDNm1sL5NGjR2nQoAHPPvssd999N7t37y7xz2TPYlEHOJnreapt3Q2EEGOFEEeAj4FnbasbAVIIsVYIsUMI8Yodc5Ypj/VogFYj+GrDEUdHURQFGDJkCHq9nqSkJAICApgzZw6jR4/m6NGjtGjRgoceeogFCxbk2YqYN28eY8eORa/X33Dq6Omnn2bBggV07tyZQ4cO3TBRkr+/P02bNmXUqFHX1n333Xe0aNGCNm3acPDgwRvuzCop4ubzaCV2YCEGAX2klGNsz0cAHaWUz9xi+6G27R8WQrwEjAU6AFlANPA/KWX0Tfs8DjwO4O/vH7p06dIi583MzMTb27vI+9tb7nwL9+ew4aSJyWEeVHF3nk74Zel36KxUxtvj5+dHSEjIDetK8gK3vRQnY1ZWFp07d2bTpk34+fnd1r7Jycmkp6ffsC4iIiJBStm+wJ0LM0NSURZAD6zN9XwCMCGf7TVAuu3xQ8D8XK+9Abyc3/tVpJnyTp6/LIMnrJJvrtjruEB5KEu/Q2elMt6e/fv3/2fdpUuXHJDk9hQ147p162TdunXl1KlTi7R/Xr8vnGCmvG1AQyFEkBBCZysAK3NvIITIfSNwf+Cw7fFaoJUQwtN2sTsMyH1hvEILqOzJfe3qsOSvE/ybkePoOIqilJI77riDEydO8Pzzz5f6e9utWEgpTcA4rF/8B4BlUsp9QohJtjufAMbZbo1NBMYDD9v2vYD1zqhtQCKwQ0q5yl5Zy6KnwkMwmi3M3nzU0VEURakA7NopT0q5Glh907qJuR4/l8++i7DePqvkIaiaFwNa12ZR/HGe7BFM5eJ2/FMURcmH81wdVW7b2IgQLhvMzItLcXQURVHKOVUsyrBG/j70bV6T+X8e41J22emzqChK2aOKRRk3rmcIl7JNfBt/3NFRFKXCql+/Pi1btqRNmza0b2+9C/X777+nefPmaDQatm/fXqjjpKSk0KJFC3tGLTJVLMq4FnX8iGhcnTmbj5FlyHtODEVR7C8mJobExMRrhaFFixb8+OOP9OjRw8HJSoYqFuXAuJ4NOX/ZwOKtJxwdRVEUm6ZNm9K4ceMCt0tISKB169bo9XpmzJhxbX1KSgrdu3enXbt2tGvXjri4OMA6NtSKFSuubTds2DBWrlzJvn376NixI23atKFVq1YcPnz4P+9VHGqI8nIgNLAyXYKr8s3GowzvHIi7q3P3XlUUu1jzGpzZg4fZBNoS+mqr2RL6fVjgZkIIevfujRCCJ554gscff7zQb1FWhilXLYtyYlzPEM5m5PB9QqqjoyhKhfPnn3+yY8cO1qxZw4wZM9i4cWOh9itLw5SrlkU5oW9QldDAysyMPcJDHeriWpRZ+RSlLLO1AK44YIjy2rVrA9bWwL333stff/11y2sVo0aNIiEhgYCAABYvXlxmhilX3yjlhBCCcT1DOHXxCj/tPOXoOIpSYVy+fJmMjIxrj3///fd872iaN28ef/75J6tXry5Tw5SrYlGOhDeqTos6vnwZk4zJbHF0HEWpEP755x+6detG69at6dixI/3796dv37789NNPBAQEEB8fT//+/enTp0+e+5eZYcoLM9pgWVgq0qiz+Vmz528Z+Oqv8uedqfYNlIfy8jt0JJXx9lS0UWevunz5smzQoIG8ePHibe3nrKPOKg7Qu5k/jfy9mb4+GYvFPnOVKIriOH/88QdNmjThmWeeue35LIpDXeAuZzQawdiIEJ5bmsjv+8/Qt0UtR0dSFKUEXR2mvLSplkU5dFer2gRV8+KL9clXJ49SFEUpFlUsyiGtRvBUeDD7Tl8iNulfR8dRFKUcUMWinLq3bR3qVPLg8/WHVetCUZRiU8WinHLVangyPJidJy4SfyTN0XEURSnjVLEoxwaFBlDDx40v1ic7OoqilFtJSUm0adPm2uLr68u0adPYtWsXer2eli1bMmDAAC5dulTgsdQQ5YpDuLtqebxHA+KPprE95byj4yhKudS4cWMSExNJTEwkISEBT09P7r33XsaMGcOHH37Inj17uPfee5k8ebKjoxaLKhbl3NBO9ajipWN6jGpdKIq9RUdHExwcTGBgIElJSdfGh+rVqxfLly/Pcx81RLniFDx1LjzaLYjJa5PYk5pOy4DS68SjKKXpo78+4uD5g5jNZrTakhmmv0mVJrza8dVCb7906VKGDBkCWCc/WrlyJQMHDuT777/n5MmTee6jhihXnMZIfSC+7i5MjynZvzQURbnOYDCwcuVKBg0aBMDcuXOZMWMGoaGhZGRkoNPp/rOPGqLcRgjRF/gM0AKzpZQf3vT6k8BYwAxkAo9LKffner0esB94S0r5iT2zlmc+7q480jWIz6MPk3Qmg8Y1S3f4ZkUpDVdbABkOGKIcYM2aNbRr1w5/f38AmjRpwu+//w7AoUOHWLVqFaCGKP8PIYQWmAH0A5oBQ4QQzW7abLGUsqWUsg3wMfDpTa9PBdbYK2NFMqpLfbx0WmaoaxeKYhdLliy5dgoK4OzZswBYLBbeffddnnzySUANUZ6XjkCylPKolNIALAUG5t5ASpn7XjIv4FrvMSHEPcBRYJ8dM1YYlb10DNcH8uvu0xz9N9PRcRSlXMnKymLdunXcd99919YtWbKERo0a0aRJE2rXrn3DcOK5lZUhyu15GqoOkPuKTirQ6eaNhBBjgfGADuhpW+cFvAr0Al6yY8YKZUy3Bsz/M4WvYo8weVBrR8dRlHLD09OTtLQbO78+99xzPPfccwXuGxoayq5du649f+uttwBo2LDhDa2DDz744NrjrKysaxe9r5owYQITJkwo6kcokD2LRV4n4v4z7oSUcgYwQwgxFPgf8DDwNjBVSpl5q/N5AEKIx4HHwVppY2Njixw2MzOzWPvbW0nl61FHw487UunolUZ1z5JtWFaU36E9qYy3x8/P79osdVeZzeb/rHM2xckYExPD2LFjGTt2LBqN5raOk52dXfT/doWZ9KIoC6AH1uZ6PgGYkM/2GiDd9ngTkGJbLgLngXH5vZ+a/KhwTl/MkiH/t0rO/OBbKd9/X8q4uBI5rpQV53doTyrj7amokx8VVXEmP7Jny2Ib0FAIEQScAh4ChubeQAjRUEp59X7O/sBhACll91zbvAVkSimn2zFrhVHLz4Px3ucZ+cajSGlG6HQQHQ16vaOjKUqRSClveUeRcp0s5oCidrvALaU0AeOAtcABYJmUcp8QYpIQ4m7bZuOEEPuEEIlYr1s8bK88ynVDso7gajYhzGYwGMBJTikoyu1yd3cnLS1NjaxcACklaWlpN9x+e7vs2s9CSrkaWH3Tuom5Hhd49UdK+VbJJ6vYKt3ZG8PHH2AyGtHqdIjwcEdHUpQiCQgIIDU1lX//vT5vS3Z2drG+FEuDIzK6u7sTEBBQ5P3VcB8VkV7PkSUr+OWzxbQcNpB+6hSUUka5uroSFBR0w7rY2Fjatm3roESFUxYy3kwViwqqyb29ee20F2sumehjkWg06pyvoii3psaGqqCEEDzavQHHzl1m/cGzjo6jKIqTU8WiAuvXoia1/NyZvfmoo6MoiuLkVLGowFy1Gh7pUp8tR8+z91S6o+MoiuLEVLGo4B7qWA9PnZa5m485OoqiKE5MFYsKzs/DlcHt6/LL7tP8cynb0XEURXFSqlgojOpaH5NFsjA+xdFRFEVxUqpYKARW9aJXU3+itp7giqFkp2JUFKV8UMVCAeDRbkFczDKyfEeqo6MoiuKEVLFQAOgYVIWWdfyY++cxLBY1zo6iKDdSxUIBbJ30ugVx9N/LxB5SnfQURbmRKhbKNXe2rEVNX3fmqNtoFUW5iSoWyjU6Fw0juwTyZ3Ia+09fKngHRVEqDFUslBsM7VgPD1ctc/9UrQtFUa4rsFgIIbYLIcYKISqXRiDFsSp56nggNICViac5m6E66SmKYlWYlsVDQG1gmxBiqRCij1BzGJZro7rWx2ixsCj+uKOjKIriJAosFlLKZCnl60AjYDEwFzghhHhbCFHF3gGV0tegujeRTWqwaOsJso2qk56iKIW8ZiGEaAVMASYDy4EHgEvAevtFUxzp0W4NOH/ZwE87Tzk6iqIoTqDAmfKEEAnARWAO8JqUMsf20lYhRFd7hlMcp3ODKjSr5cuczcd4qENd1JlHRanYCtOyGCSljJRSLr5aKIQQQQBSyvvsmk5xGCEEY7oHkXw2kw2H/nV0HEVRHKwwxeKHQq5Typm7WtWmho+b6qSnKMqtT0MJIZoAzQE/IUTuFoQv4G7vYIrj6Vw0PNylPpPXJpF0JoPGNX0cHUlRFAfJr2XRGLgLqAQMyLW0Ax4rzMGFEH2FEElCiGQhxGt5vP6kEGKPECJRCLFZCNHMtr6XECLB9lqCEKLn7X4wpWQM7VgPd1eNmklPUSq4W7YspJQrgBVCCL2UMv52DyyE0AIzgF5AKtZ+GiullPtzbbZYSjnTtv3dwKdAX+AcMEBKeVoI0QJYC9S53QxK8VX20nF/uwC+T0jl5b6Nqebt5uhIiqI4wC1bFkKIV2wPhwohPr95KcSxOwLJUsqjUkoDsBQYmHsDKWXuAYi8AGlbv1NKedq2fh/gLoRQ31IOMrpbEAaThUVbVCc9Ramo8rt19oDt5/YiHrsOcDLX81Sg080bCSHGAuMBHZDX6ab7gZ25btnNve/jwOMA/v7+xMbGFjEqZGZmFmt/e3N0vtbVtczZeJhm4hQ6bd630To6Y0GcPR+ojCXB2fNB2cj4H1JKuyzAIGB2rucjgC/y2X4osOCmdc2BI0BwQe8XGhoqiyMmJqZY+9ubo/NtPvyvDHz1V/ndXyduuY2jMxbE2fNJqTKWBGfPJ6VzZQS2y0J8p+d3N9Qv2E4L3aLI3F1AHUoF6uZ6HgCcvsW2YD1N9VWu9w8AfgJGSimPFPBeip11Ca5Kk5o+zN58lEHtA1QnPUWpYPI7DfVJMY+9DWho68B3CuuAhENzbyCEaCilPGx72h84bFtfCVgFTJBS/lnMHEoJuDqT3ss/7GZz8jm6N6zu6EiKopSi/O6G2lCcA0spTUKIcVjvZNICc6WU+4QQk7A2e1YC44QQdwBG4ALwsG33cUAI8IYQ4g3but5SSjXfpwPd3aY2H/2WxOxNx1SxUJQKJr/TUMuklIOFEHvI43SUlLJVQQeXUq4GVt+0bmKux8/dYr93gXcLOr5SutxctIzUB/LpukMc/ieDhv6qk56iVBT5dcq7+kV+Fzd2yru6KBXQsE71cHPRqJn0FKWCuWWxkFL+bXv4tJTyeO4FeLp04inOpqq3G/e1q8OPO06Rlvmfu5kVRSmnCjOQYK881vUr6SBK2TG6axA5JgtRW084OoqiKKUkvx7cT9muVzQWQuzOtRwDdpdeRMXZNPT3IaxRdRbGHyfHpGbSU5SKIL+WxWKs1yZWcuO1ilAp5fBSyKY4sTHdgziXmcPKxPy6ziiKUl7kd80iXUqZIqUcYrtOcQXrXVHeQoh6pZZQcUrdQqrR2N+HOZuPXe1tryhKOVbgNQshxAAhxGHgGLABSAHW2DmX4uSudtI7eCaDuCNpjo6jKIqdFeYC97tAZ+CQlDIIiARUr2qFu9vUppq3Ts2kpygVQGGKhVFKmQZohBAaKWUM0MbOuZQywN1Vy/DOgaw/eJbks5mOjqMoih0VplhcFEJ4AxuBKCHEZ4DJvrGUsmJ450B0LhrmqU56ilKu5TeQ4FUDgWzgBWAY4AdMsmcopeyo5u3GvW3qcGTlOoaLA+DmBnq9o2MpilLCCiwWUsrLuZ4usGMWpYx62v0sNRZNQGc2wdIoiI5WBUMpH+LjITYWwsMr/P/T+Q0kmMGNAwgK23MBSCmlr52zKWVE4J5tmC0mtNKCNBgQsbEV/h+WUg7Ex2PuGQk5OWjc3RAV/I+g/PpZ+EgpfXMtPrl/lmZIxcmFh6Nxc8MkNORoXMjp1t3RiRSlWM5l5rDq8yXInBy00gIGg7WFUYEV5poFQojWwNVvgI1SSjXch3KdXo+IjuavL2bzibYpbdOr8EbBeymK0zFbJIu3Hmfy2iQa6wLprdMhTUaETmc9FVWBFVgshBDPAY8BP9pWRQkhvpFSfmHXZKUoPj6eqKgo3Nzc0FfgZmax6PUYcnJokV6NOZuP0bNJDbqGVHN0KkUptJ0nLvDGir3sPXWJriFVefvpJ3B9sou6ZmFTmJbFo0Cnqxe6hRAfAfFAuSgWcXFxRERGYDKYiIqKIjo6WhWMYpjQrymbD5/jpe938dvzPfDzcHV0JEXJ14XLBj5ee5Cl205Sw8eNL4a05a5WtazzzNfQV/gicVVh+lkIIPfQombbunLhp7U/YTAYsFgsGAwGYiv4ecni8tBpmfpgG85m5PDmir2OjqMot2SxSJb+dYKeU2JZtj2VMd2CiH4xnAGta1sLhXKDwhSLecBWIcRbQoi3gC3AHLumKkX39b0PnU4HGhAugrCwMEdHKvNa163Esz0b8nPiaX7drUalVZzP3lPp3PdVHK/9uIeGNXxY/Wx3Xu/fDG+3Ql3GrZAK08/iUyFELNANa4tilJRyp72DlRa9Xk/s+lhe/fJVTtU/xQGfA3Shi6NjlXljI4JZn3SW13/aS/vAKtT0c3d0JEUh/YqRKb8nsWjLcap46fh0cGvubVtHtSQKoTAXuD8DvpNSfl4KeRxCr9fzdvbbrNGsYdqOadTzrUevwLwmCFQKy0WrYerg1tz5+SZe/mEXC0d3VP8gFYeRUvLjjlN8sOYA5y8bGKmvzwu9GqlrarehMKehdgD/E0IkCyEmCyHa2zuUIwgheKfbO7Su3pr/2/R/7D2nzrcXV4Pq3rzevxmbDp/j2y3HHR1HqaAOnrnE4K/jefH7XdSt4snKcd146+7mqlDcpgKLhZRygZTyTqAjcAj4yDa/RYGEEH2FEEm2QvNaHq8/KYTYI4RIFEJsFkI0y/XaBNt+SUKIPrfxmYrMTevGZxGfUdWjKs+sf4a/M/8ujbct14Z3qkdYo+q8v/qAGplWKVUZ2Ube+XU//T/fTPLZTD66vyXLn+xCizp+jo5WJhWmZXFVCNAEqA8cLGhjIYQWmAH0A5oBQ3IXA5vFUsqWUso2wMfAp7Z9mwEPAc2BvsCXtuPZXVWPqkzvOZ1sUzZj14/lsvFywTsptySEYPIDrXB31TJ+WSJGs8XRkZTyLD6eulFRbFq4ksgpG5j75zEe7FCX9S+G82CHemg06lRoURVmpryrLYlJwF6sc3APKMSxOwLJUsqjUkoDsBTrCLbXSCkv5XrqxfWxqAYCS6WUOVLKY0Cy7XilIqRyCFPCpnD04lFe3vAyJosakb04avi688G9Ldmdms709cmOjqOUV/HxWHpGUn/OXNo/Oogeacn89HRX3r+3JZW9dI5OV+YVpmVxDNBLKftKKedJKS8W8th1gJO5nqfa1t1ACDFWCHEEa8vi2dvZ15661OnC/3X6Pzad2sQn2z8pzbcul/q1rMV97eowPSaZnScuODqOUo5cMZhZtu0kCz+Yj8U2lpObxczH1c7Tpm4lR8crNwpz6+zMIh47r/ae/M8KKWcAM4QQQ4H/AQ8Xdl8hxOPA4wD+/v7F6lCXmZn5n/1rUIMInwiiDkRh+MdAD58eRT5+ceWVz9kUlPGOypINOnhyfjyTunjg5lK6pwTKw+/QGThLxlOZFmJPGtl8ysQVE/Su2YiHXF2xmIxIVxd2+flxyQly5sVZfoe3RUpplwXQA2tzPZ8ATMhnew2Qnte2wFqsrZtbvl9oaKgsjpiYmDzXm8wmOe6PcbLVglZyU+qmYr1HcdwqnzMpTMa45HOy/mu/ytd/2m3/QDcpL79DR3NkxmyjSa5IPCUHz4yTga/+KkP+b5V8ZvEOueXIOWmxWKSMi5NHxoyRMi7OYRkLw5n+OwPbZSG+0+3ZXXEb0FAIEQScwnrBemjuDYQQDaWUV++s6g9cfbwSWCyE+BSoDTQE/rJj1lvSarR81OMjHv7tYV7a8BIL+y2kUeVGjohSLuiDqzKmWxCzNh0jsqk/EY1rODqSUgacPJ/F4r9OsGzbSdIuG6hbxYNX+zZhUPsAqnm7Xd9Qr+dETg4N1HhOJS6/yY+q5LejlPJ8Aa+bhBDjsLYKtMBcKeU+IcQkrJVsJTBOCHEHYAQuYD0FhW27ZcB+rPN9j5VSmvN8o1Lg6erJFz2/YOiqoYyLHsfi/oup5qFGVC2qF3s3ZuOhc7zyw27WPt+DKurio5IHs0Wy/uBZorYeZ8OhfxFAZFN/hncOpHtINXVnUynLr2WRwPWZ8eph/TIXQCXgBBBU0MGllKuB1Tetm5jr8XP57Pse8F5B71FaanrV5IvILxj12yieXf8sc/vMxd1FDWFRFO6u1sEGB87YzOs/7eHLYe1U727lmn8uZfPdtpMs/esEp9Oz8fd149meDXmoY11q+Xk4Ol6Fld9MeUFSygZYWwYDpJTVpJRVgbu4PrdFhdK8anM+6P4Be8/t5fXNr2ORqs9AUTWr7cuLvRuzZu8Zftp5ytFxFAezWCSbD5/jqUUJdPlwPZ+uO0RwDW9mDg9l86s9eaFXI1UoHKww1yw6SCmfvPpESrlGCPGOHTM5tch6kYwPHc+UhCkE7gzk2XbPFryTkqfHujdg/YGzvLliHx2DqhBQ2dPRkZRSdmn9Rg4u/YUFuvqs8q5PZU9XHu0WxNCO9ahfzcvR8ZRcClMszgkh/gcswnpaajiQZtdUTu7h5g+TcimFWXtmEegbyMCQgQXvVAwbN29k3rfzyt1MflqNYMrg1vT7bBMvLtvFksc6q/PQFUTSmQyi5/7MqDfH0M5sorWLK4O/WUanYX1xdy2VwRqU21SYTnlDgOrAT7alum1dhSWE4PXOr9OpZifein+LbWe2lfh7ZJuyiT4RzcivRhIRGcHC+QvtJztFAAAgAElEQVTpGdmT+Pj4En8vR6pbxZOJA5qx9dh55mw+5ug4ih1ZLJI/9v/D0Flb6DNtI1fWRaOzmHCRFtwsJsL+3qcKhRMrTKe888BzQghvKaUaCc7GVePKlPApjFgzghdiXyDqzigCfQOLdcwsYxabT21m3fF1bEjdwBXTFTI2ZCBNEiyQY8ghNja2XLUuAAaFBvDH/n+YvDaJ7o2q0aSmr6MjKSUoI9vI99tTWRCfwvG0LGr5ufNK38YM7z0a7V3LwGAAnc46z7XitAozn0UXYDbgDdQTQrQGnpBSPm3vcM7Oz82PGT1nMHT1UMZGjyXqzij83G5vRMvLxstsTN3IuuPr2JS6iWxzNlXcq3BXg7voFdgLU4iJPiv7kJ2TDVpo0bGFnT6N4wgh+OC+lvSZtpHnlyayYlxX3FzUX5hl3fG0y8yPS+H77alk5phoV68SL/dpTJ/mNXHVaoAQiI6G2FhroShnfwSVN4W5ZjEV6IO1oxxSyl1CCMeNe+Fk6vrW5bOIzxjz+xiej3meb3p9g6s2/3HyMwwZxJ6MZd3xdfx56k8MFgPVPKpxT8g99K7fm3Y12qHV2L4sa0N0dDQzZs9gW51tJHgkMIDCjONYtlT1duOj+1vx6ILtTF13mNf6NXF0JKUIpJTEH0lj7p/HiD54Fq0Q3NWqFqO6BtE6r3Ga9HpVJMqIQvXgllKevOk+eId1kHNG7fzbManrJCZsmsCkLZOY1GXSf/oNpOekE3MyhnXH1xF3Og6TxYS/pz+DGw+mV2Av2tRog0bkfQlJr9eTk5NDM49mLEtaxugWo6nnW680Plqpimzqz5CO9fh64xF6NqlBx6B8+4UqTiTbaObnnaeY92cKSf9kUMVLx7iIEIZ3DsTfV/VHKg8KUyxO2k5FSSGEDuvIsAfsG6vsuavBXRy/dJyZu2ZiOmbC96Qv7fTtyKidwR/H/2Dr31sxSRO1vWozrMkwetXvRctqLW9ZIPLyWMvH+OnwT8zcNZP3u79vx0/jOP/r35S4I+eYMzmK1rXScbsjUv3l6cTOpGfz7ZYUFm89wYUsI01r+fLxA624u3VtdbG6nClMsXgS+AzrEOGpwO/AWHuGKquebv00cXFxfPTCR0iTRLgIgl4JonHbxoxsPpLegb1pVrVZkXsrV/eszpAmQ5i/bz6PtnyU4ErBJfwJHM/LzYWvgw0E/t9LuFhMyA/eR0RHq4LhLOLjqRcVRdKFHKZn12DNnr8xS0mvpv6M6hpE5wZVVG/8cqowd0OdA4aVQpYyTwhBkwtNrKNZWUCYBQMZyJR7p5TYP6BRLUbxXdJ3fJn4JVPCp5TIMZ1Nk6QdWCwmNBYL5hwDmpgYhCoWDmf5Mw4ZGUmgwYBh/kIujPiQh+/rw8P6+tSrqjpUlne3PAcihHjF9vMLIcTnNy+lF7FsuaPnHbi7uaPVanHTuTHozkEl+pdWZffKjGg2gt+P/86BtHJ6NjA8HOHmhlmjxaDRMscl8OpQ9YqD7Dp5kUUfLUAaDNcmF5oTlMUbdzVThaKCyK9lcfWbaHtpBCkv9Ho90dHRxMbGEh4ebpc+ESObj2TxwcXMSJzB9MjpJX58h9PrEdHRaGJiWKKrz7vn/Dj1634m3lX0U3hK0Zy/bGDy2oMs3XaSiFrNGarTYTEa0LjpcLujp6PjKaXolsVCSvmL7eeC0otTPuj1ert2nPPV+TKq+Sg+3/k5u/7dRevqre32Xg6j1yP0ekZJyclf9zPvzxQ0QvC//k1VwSgFZotkyV8n+OT3JDKyTYzuGsTzd/TG5ZEOHJ07lwajR6vrSBVMfvNZ/EIeU5leJaW82y6JlEIZ1nQYiw4sYvrO6czqPcvRcexGCMHEu5ohJczZfAyNgP+7UxUMe9px4gITV+xl76lLdAqqwqSBLWhc08f6oppcqMLK7zTUJ6WWQrltnq6ePNriUSZvn8y2M9voULODoyPZjRCCNwc0wyIlszYdQyMEr/VrogpGCTuXmcNHaw7yfUIq/r5ufD6kLQNa1VK/ZwXI/zTUhquPbf0rmmBtaSRJKQ2lkE0pwODGg1mwbwHTd05nft/55foftRCCt+9ujkVKvt54FCEEr/ZtXK4/c2kxmS1EbT3BlN+TyDKYeaJHA56JbIi3mz1nXVbKmsKMDdUfmAkcwTpTXpAQ4gkp5Rp7h1Py5+7izuOtHufdre8SdzqOrnW6OjqSXQkhmHR3C6SEmRuOoBHwch9VMIpjW8p5Jq7Yx4G/L9EtpBpv3d2ckBrejo6lOKHC/OkwBYiQUiYDCCGCgVWAKhZO4L6G9zF371y+2PkFXWp3KfdfnBqN4J2BLbBI+DL2CBoheLF3o3L/uUva2YxsPlx9kB93nqK2nztfDmtHvxY11e9RuaXCFIuzVwuFzVHgrJ3yKLfJVevKk62fZGLcRNafXE9kvUhHR7I7jUbw3j0tkFIyPSYZjUYwvlcjR8cqE4xmCwvjjzNt3SGyTWaeDg9mXM8QPHXqlJOSv8L8H7JPCLEaWIb1msUgYJsQ4j4AKWWFnI/bmQwIHsCcvXOYvnM6EXUjbmu8qbJKoxG8f29LLFLyefRhNAKev0MVjPxsOZrGmyv2kfRPBmGNqvPmgGY0qK5OOSmFU5hi4Q78A4TZnv8LVAEGYC0eqlg4mIvGhadbP82rm15lbcpa+gX1c3SkUqHRCD68rxUWCdP+OIxA8NwdDR0dy+n8cymb91YdYOWu09Sp5MHXI0Lp3cxfnXJSbkthxoYaVdSDCyH6Yh2EUAvMllJ+eNPr44ExWEdT+hcYLaU8bnvtY6A/1iFJ1gHPSTXmwy31DerLrD2z+DLxS3oF9sJFUzFOK2g0go/ub4WUMPWPQ2gEPBOpCoaUkovRGzm87BemmeqwvXYTno1syFNhwXjo1Giwyu0rzN1QQcAzQP3c2xfUKU8IoQVmAL2wjla7TQixUkq5P9dmO4H2UsosIcRTwMfAg7Yh0bsCrWzbbcbasokt3MeqeDRCw7i243g+5nl+Pfor94Tc4+hIpUarEXz8QCuklExZdwiNRjA2IsTRsUpFZo6JlHOXOXruMkf/zeTYucscO3cZ3x3bmPXta7Qzm5jn6sr5Fauppa7rKMVQmD8/fwbmAL8Alts4dkcgWUp5FEAIsRQYCFwrFlLKmFzbbwGGX30J6+kvHdbbdV2xngpT8tGzbk+aVW3GzF0z6R/Uv8AZ+8oTrUYweVBrLFIyeW0SQsDT4eWjYBjNFk6cz+LYv9ZCkLswnM3IubadEFDbz4MG1b142HgcN4sJjbTgYjZRa+dW6KvGclKKrjDFIltKWZRRZusAJ3M9TwU65bP9o9hux5VSxgshYoC/sRaL6VLKcjrEaskRQvBM22d46o+n+Cn5JwY3HuzoSKVKqxFMGdwGCXz8WxIaIXgyrOzM+ZG1YRPus75lbfJ5/qrZ5For4cT5LMyW62dgq3jpCKrmRY9G1WlQ3YsG1bwIquZNYFXP6xMONbPA8m/AYACdzjrHtaIUgyjoMoAQYijQEOukR9f+jJFS7ihgv0FAHynlGNvzEUBHKeUzeWw7HBgHhEkpc4QQIVivdTxo22Qd8KqUcuNN+z0OPA7g7+8funTp0nw/S34yMzPx9nbeO0MKm09KybR/ppFmSmNi7YnoNLpSSGflLL9Ds0Xyze4ctp4x82BjHf2CrC0sZ8mXF8P2PfSY8BIuZhNGrQsPD3mPvxs3p6aXoKanhppeAn8vDTU9NXjrCndh2nffPiolJnKxTRsuNW9eYlmd+fcIzp8PnCtjREREgpSyfYEbSinzXYAPsLYKNgAxtmV9IfbTA2tzPZ8ATMhjuzuwDodeI9e6l4E3cj2fCLyS3/uFhobK4oiJiSnW/vZ2O/n++vsv2WJ+C7lg7wL7BcqDM/0OjSazHBuVIANf/VXO2nhESulc+XJLOH5efhY5ShqFRkqQFq1Wmt97z9GxbslZf49XOXs+KZ0rI7BdFvB9LqW89eRHudwLNJBShkkpI2xLYU5+bgMaCiGCbGNLPQSszL2BEKIt8DVwt5Qyd0e/E0CYEMJFCOGK9eK2Og1VSB1qdqBTrU7M2TuHLGOWo+M4hItWw7QH29C/ZS3eXXWAOZuPOTpSnmIOnmXorC0cbNIOjZsbFo0GodOhiYhwdDRFuUFhisUuoNLtHlhKacJ6amkt1i/6ZVLKfUKISUKIq3dSTQa8ge+FEIlCiKvF5AesY1Htsb3/LmmbX0MpnGfaPsP57PMsPrjY0VEcxkWrYdpDbbizZU1Wfb2cnOkLIT7e0bGuWZ6QypiF2wmp4c3bHzyGZn00KaNHg5pzXHFChbnA7Q8cFEJs48ZrFgXOZyGlXA2svmndxFyP77jFfmbgiUJkU26hdfXWhAWEMXfvXAY3HoyvztfRkRzCVavh88BsLMv+h8ZoxPDLYizr1uHeo7vDMkkp+WbjUT5Yc5CuIVWZOTwUH3dXNVeE4tQK07J4E+upqPexDip4dVGc3Ng2Y8kwZPDt/m8dHcWhXDZtxNVswkVa0BiNfPvhQnanXnRIFotF8u6qA3yw5iD9W9Vi7iMdrIVCUZxcgcVCSrkhr6U0winF07RqU3oF9uLb/d9yIfuCo+M4Tng4QqezXg9w07G1Xivu+zKOGTHJN9ySam8Gk4UXliUyZ/MxHulSny8eaoubi+pNrZQNtywWQogMIcSlPJYMIcSl0gypFN3YNmPJMmYxb988R0dxHL0eoq3XA7Tr1/PJp0/Qp3lNJq9NYsisLZy6eMXuES7nmHh0wTZWJJ7m5T6NeXNAMzQaNTaTUnbcslhIKX2klL55LD5Syop5ArwMCq4UTP8G/VlyYAnnrpxzdBzH0es5MWwY6PVU8tQxfWhbPhnUmn2n0uk7bSO/7Dptt7dOy8xh6Kwt/Jl8jo/ub8nYiBA1iJ9S5pT/sawVnmr9FEaLkdl7Zjs6itMQQvBAaACrn+tOcHVvnlmyk/HLEsnINpbo+5w8n8UDM+M5eCaDr0e058EO9Ur0+IpSWlSxqADq+dbjnpB7WJa0jL8z/3Z0HKcSWNWL75/U82xkQ37eeYo7P99EwvHzJXLsA39f4v6v4kjLzCFqTCd6NfMvkeMqiiOoYlFBPNHKeify17u/dnAS5+Oq1TC+VyOWPaFHShg0M56p6w5hMt/OuJk32nI0jcEz49EIwQ9PdaF9/SolmFhRSp8qFhVELe9aDGo0iJ+Tf+bEpRMlfvz4+HiioqKId6JOb7erff0qrH6uO/e0qcNn0YcZ/HU8J9Juvwf8b3v/ZuTcv/D3c2f5011o5O9jh7SKUrpUsahAxrQcg6vGlZm7Zhb7WFJK0q6ksfvf3Uz5YQphPcOYO3cukZGRZbpg+Lq78umDbfjsoTYcPptJv8828kNC6tUxygq0aMtxno7aQfPavnz/hJ46lTzsnFhRSkfFmE5NAaC6Z3WGNBnC/H3zebTlowRXyn/47ixjFqcyT11bUjNSSc1MJTUjlVOZp7hist5y+u+v/2LMMYKEHEMOsbGx6Mt4L+SBbeoQGliZ8ct28dL3u4hJOsv797TEzzPvDnRSSj6LPsy0Pw7Ts0kNpg9ti6dO/fNSyg/1f3MFM6rFKL5L+o7/Rf2PFuktaNmpJf7N/G8oBlcfn8++8UKvp4sndXzqUNenLp1rdSbAJ4AA7wDO1jzLyFUjyc7JRmol52ufR0pZ5m8PDajsyZLHOjNzwxGmrjvEjuMX+HRwG/TBVW/YzmyRTFyxl6itJ7i/XQAf3t8SV61qtCvliyoWFUxl98p0M3Zj6vipLDUtRbgIgl4JwjPEE63QUtOrJgE+AUTUjSDAJ4A63nUI8A6gjk8dKrtVzrsA1IWA6ABmzZ2FOdTMb/yG2CR4p+s7uGndSv9DliCtbYrWbiHVeP67RIbO3sITPYIZ36sROhcN2UYzzy9N5Ld9Z3gyLJhX+zYu80VSUfKiikUF5JPqAybAAsIs6GPuw+v3v46/pz8umqL9L6HX68nJySEsLIw5e+fw2Y7POJ15mmkR06jmUa1kP4ADtK5biV+f6cY7v+5n5oYjbE7+l0/rXiFhwU+c9Q7mjSfu59FuQY6OqSh2o9rKFVCfyD64u7mj1Wpx07kxYsAI6njXKXKhyE0IwZiWY/g0/FOSzicxbNUwDl84XAKpHc/LzYUP72/FzOGhVNmdQN1BAxi04huWff8Gj2rPODqeotiVKhYVkF6vJzo6mnfeeYfo6Gi7XIzuFdiL+f3mY7QYGbFmBBtTNxa8Uy7xJ+P5YNMHxJ90vjur+raoyYzaGegs1pFsXUxGiI11dCxFsStVLCoovV7PhAkT7HrXUvOqzVncfzH1fOrxzPpniDoQVahbUONPxhO5IJw31r9O5MKeTlkwfPregdbNDbRa0OkgPNzRkRTFrlSxUOyqpldN5vedT3hAOB/+9SHvbX0PoyX/8ZdiU2IxmI2YkRiM2cR+Pww2fgJnD0Ih+zvYnW0kW955R81sp1QI6gK3Yneerp5MjZjKtB3TmLd3HicuneCT8E9uOXtfeP1wdC7uGMwGdBoN4W6VYf071qVKMDTpD03ugoAOoHHg3zt6vSoSSoWhioVSKjRCw/jQ8QT5BjFpyySGrx7OjJ4zqOtb9z/b6uvqiR4ZTWxKLOH1w9HX1cOl05C0Gg6ugi1fQtzn4FUDGvezFo6gHuDq7oBPpigVgyoWSqm6t+G9BPgE8ELsCwxdPZRpEdMI9Q/9z3b6unprkbjKtzZ0GGNdrlyE5D/g4K+wdznsWAA6bwi5w1o4GvYCj0ql+KkUpfxT1yyUUtehZgei7oyiklslxvw+hpVHVt7eATwqQcsHYNB8eOUoDPvB+vx4HPw4BiYHw8J74K9ZkH7KLp9BUSoa1bJQHCLQN5BFdy7ixdgXeX3z6xxLP8YzbZ9BI27z7xcXN2tLomEv6D8VTm23tjgOroLVL1mX2u2Ir9WcZf/8g1sDHfp6XezzoRSlHFPFQnEYPzc/vur1Fe9vfZ/Ze2aTkp7C+93fx8OliCO1ajRQt6N16TUJ/j0EB38lfvdiIhNmYADmz1tOdPBA9K2GQFAY+KgJiRSlMOx6GkoI0VcIkSSESBZCvJbH6+OFEPuFELuFENFCiMBcr9UTQvwuhDhg26a+PbMqjuGqcWVi54m83P5lok9E88hvj3A262zJHLx6I+g+nnUt7sZ8RgdxWnJOSmKPx8KPj8GURvBlF/jt/+DwOjBcLpn3VZRyyG4tCyGEFpgB9AJSgW1CiJVSyv25NtsJtJdSZgkhngI+Bh60vbYQeE9KuU4I4Q0UfdoyxakJIRjZfCSBvoG8svEVhqwawhifMaTsTCE8PPw/HQellGQaM0m7kkZadlq+P89nnyftYBrGeUakUSJcBO5Rb0DXCDgaC0diYNts2DIDNK5QtxMEh0ODCKjdFjRah/xOFMXZ2PM0VEcgWUp5FEAIsRQYCFwrFlLKmFzbbwGG27ZtBrhIKdfZtsu0Y07FSYTVDWNhv4UM/2o4wyYNAxO4uLow/PPheAZ73lAEDBbDf/YXCCq7V6aKexWqelSldfXWVPWoyta/tnLMdAwkSLPks+Wf0zC0Ef27Po/o9gIYr8CJeGvhOBoL69+1Lu5+1ltyG4Rbi0eVBqBGlFUqKFHYGcBu+8BCPAD0lVKOsT0fAXSSUo67xfbTgTNSyneFEPcAYwADEAT8AbwmpTTftM/jwOMA/v7+oUuXLi1y3szMTLy9vYu8v705ez4ouYxzv53Lt/O/tbYlNVDv/no0vacpPloffLQ++Gp98dH4XHt+dZ2Xxgut+G9LYN++fbz44osYjUZcXF3oMKED6fXTaebRjIeqPERll8o3bO9qSKfyhV3XFvecfwHIdqvB+SqtuVC5NRcrtSbxykkS0xNp49eG5n7Ni/25oWL9d7YXZ88HzpUxIiIiQUrZvsANpZR2WYBBwOxcz0cAX9xi2+FYWxZutucPAOlAA6ytn+XAo/m9X2hoqCyOmJiYYu1vb86eT8qSyxgXFyc9PDykVquVHh4eMi4urkSOOWbMGBkXFydNZpNctH+R7LCog+y4qKNccmCJNFvMee9osUh5LlnKrd9IuWSolO/XlfJNXxn3pqf0eEsjtW8J6fGOm4xL2VjsjFJWrP/O9uLs+aR0rozAdlmI73R7noZKBXJ3zw0ATt+8kRDiDuB1IExKmZNr353y+imsn4HOwBw75lWcxNVRcWNjY/O8ZlHUY+bk5Fw71rCmwwivG87bcW/z3tb3WHNsDW91eYsgv5vmpBACqgZbl46PgdkEfycSu/4NDMd+wwwYTDnEfjsAfaOB0LC3tXOgustKKWfsWSy2AQ2FEEHAKeAhYGjuDYQQbYGvsZ6uOnvTvpWFENWllP8CPYHtdsyqOBm9Xm/3ebzreNfh615fs/LISj7e9jEPrHyAp9o8xcPNH8ZVk/dc22hdIKA94RET0Z3cYB2/SqslvEEknNgC+3+2blerta1w9IKA9upCuVLm2a1YSClNQohxwFpAC8yVUu4TQkzC2uxZCUwGvIHvbVNRnpBS3i2lNAshXgKihfWFBGCWvbIqFZcQgoEhA+lapyvvb32fz3Z8xtqUtbzd5W2aVW12y/3yHL9KSjizBw7/bh2OZNMU2DgZPCpDcE9r8QiOBO/qpfgJFaVk2LVTnpRyNbD6pnUTcz2+I5991wGt7JdOUa6r5lGNT8M/Jfp4NO9ufZehq4bycPOHear1U7i75D1A4X/GrxICarWyLj1egisXrHdYHV5nLR57lwPCektuw17W4qFuz1XKCNWDW1FyiQyMpH3N9nya8Clz984l+kQ0b+nfon3Ngm8W+Q+PytDiPutiscCZXXD4D2vLY+Nk2PAReFSxXuNo2It4D1+iTizH7aTbjUVIUZyAKhaKchM/Nz/e7vI2/YL68VbcW4xaO4rBjQbzQugLeOuKeLujRmNtRdRuC2EvQ9Z5OLL+Wqsjfs9iIsnCICAqZS7RgX3R+7cCb3/wqXnjT4/Kqr+HUupUsVCUW+hcqzM/3v0jMxJnsOjAIjakbmCifiI9AnoU/+CeVawj5bZ8ACwWYtc8i2H7l9bZAaWF2DM70Z9KBEMe/VG1brbi4Z+riNS0Pb/+M/7CYWKPb7p+TUVRikEVC0XJh6erJy93eJk+9fvwZtybjI0eS7+gfrzW8TWquFcpmTfRaAhvNQxd4lxyTDnoXNwIH7oc6uohJxMy/4GMM5B5BjL+ufFn2hE4/qf1+kgu8ZiIFFkYAJ3QEh06Fn2zB6wtGzfn6AwG1vnWb7hJoASOF3UiqmRO5UlJfMoGYo9FEx7cG31g92LnK8tUsVCUQmhVvRXL7lrG7L2z+Wb3N8Sfjmegy0AyD2YSHhZO967d0RbjQrW+rp5pzacxK2oWjw177PoXnZu3dakanP8BTDm2omItIrGJ8zAc+tnWUjETu+1r9Nvmg9BA9SZQpx3UaQ91QqFGM+stwaUhJwPSkuHcYeKPRhO5ew4GabYWtPr90PsEgNbVtuisPzW5HuezPv5CMpHrJ5BjNhB1YiHRXSei9w0AY5ZtuQKGXI+Nl20/r1jXGbKuPzZmEW+8ZD01COjiJhP9cEyFbqGpYqEoheSqdeWp1k/Rq14vnp7zNK+8/grSKJnkOomgV4LwbuiNq8YVnUaHq9YVF40LrhpX66J1vfZYp9VdX29b/jnwD4ufX4zZaGbvqr20jG55e/1MXNygUj3rAoR7V0N39DdbPxAd4YN/BLSQuh1OJcDB1bBzkW1fD6jdxlo46oRa+4X41S36dRGLGdJPwrlkSDsM5w7BucPWIpHx97XNYoUBA2Zrx0ZpJvbvBPRpKWA2gMUIZqP1sdkIN470k6dYcjCQg0WAwWQgdsM76HHLtYUAnRe4eoCrp22xPfaqDpVuXB97ZhuGlD+sBddiIjYlVhULRVEKL6RyCGGGMKJN0UgpEWZBy/SWRLaMxGgxYjQbrT9vfmx7brAYyDRmXnvNZDGRtCkJk9EEFsjOyWb87PG8U/sdwgLCcNXeooNgPvLsBwLWW3bB2ifkwjFITbAWj1MJ1pkFzdOtr3vVyFU8QqF2O/CodONpnurNchWEw7l+HgFzzvUw7n5QrZF1MMZqIVC1IVRrSHjWWXRR/a4XtKun3vJiMVuLxs1FJNfP8L8T0K1+CoPZiM7FlfB7FkLdrtcLgovbbRXA8JPx6BZuvp6vfvht/3coT1SxUJQi6BnRk/fc3sNgMKDT6Xh1yKvo2xb9r874mvFE/hRJTk4OWlctxiAjL8S+QCW3StwZdCf3hNxDkypNELfxZfeffiC5CWEdRbdKA2g1yLrOZICz+2ytjx3WWQcPrbme0bcmkZlHMEgzUcfmEC090F/9ChFaqFwfqjWEkMhrBYGqDcGrWp5f0nqa5l3Q8qLR2vqj5N3nBUBfqxXR1ZswN2YuoyNGF7sVcMuCW0GpYqEoRVDS41ddPd7cuXMZPXo0HTt1JP7veFYkr+CHQz+w+OBiGlVuxMDggfRv0J+qHlVL6JPk4qK7fnvvVdnptsKRQOyebzFk2E4bAbEh4ehDn7K2GirXt+5/m/ItaEWgr6snp15OiR2zpPOVZapYKEoRlfT4VTcPdtitTje61elGek46vx37jRVHVjB5+2SmJkyle0B3BoYMpEdAj1uPY1US3P0gOAKCIwgP6opuYeT1O7bC/nfr00ZKuaOKhaI4OT83Px5s8iAPNnmQIxePsCJ5Bb8c/YWYkzFUca9y7TRV4yqN7Zrj6mmZkjrNo5QtqlgoShkSXCmY8e3H82y7Z4k7HcfPyT/zXdJ3LDqwiCZVmnBPyD3cGXQnld0rF3ywQrJIC1nGLDIMGZw+eBq5SUJDbpyAQCn3VLFQlDLIReNCj4Ae9AjowcXsi6w+tpoVR1bw4V8f8sn2TwgLCKNhekMu7L9A9x7daR7anAxDBpmGTDIMGWQYb3x882u5nzBw0hwAAA6LSURBVGcaM5FIspKzOPbxMaRJsvDbhSz5ZQn397rf0b8KpZSoYqEoZVwl90oMbTqUoU2Hcuj/27v34KjKNI/j3yfd6U4gCUGQgMglRMNFS5AoKCuiFS9cFNBlVlfAnVUGxZnRHUZHdF1nBnfYcabUgh3Z8VaL6yIyLqJhA4IggbJgWOQiiFBAQgQFVGSUCSRNd3j2j3OITUzSnUt3DuT5VHVxzum3z/lxEnj63N73L7t5Z+87vP7u68z51zloRBG/kPuLXNpd1K7Oz6dIChmpGWQGMskMZJKRmsEFGReQFcg6Y/mKbSsoj5Sjp5RwOMyPX/wxi04uYmzeWEb2Hkl2WnaS/+YmmaxYGHMOye+YzyNXPkLovRDbqrehp5znQIYcH8Kk4ZPOKAinp9v528V1S27+hHwW/8diQqEQgWCAH43/ETuqd/CbDb/h6Y1PM+LCEYzNG8vw7sOb9GyI8TYrFsacgwqvL2RWYFbNcyAP/uBBru7Tsrf3nr5ra9fRXRSVFlFcVsyq/avIDmYzKncU4/LGMaDTgEY9G2K8y4qFMeegRIxjfnq90bf3AvQ7rx/9zuvH9ILprDu4jqLSIhbtXsSCXQvo06EPt+bdyi19bqFr+64tksG0DisWxpyjkjGOebToi+7HTh5jRfkKlpQuYfbm2czZPIeh3YYyNm8shT0LaZda9/UT411WLIwxLS4rkMWE/AlMyJ/AgWMHWFK2hKLSIh7/4HHS/enc2OtGxuWNI7wvzNo1a1v06MckhhULY0xC9cjqwQODHmDawGls+XILRaVFLC9fzhvL36D8d+VoRPGl+pg8ezJ5A/MI+oIE/UHSfGkEfUHS/GkEfIGa+aA/6CyPmk/zOW22btzK/PnzCQaDVnxamBULY0xSiAiDcwYzOGcwM4bMYNpj09gX2YeeUiLhCGvXrGVPhz2EqkNURapQtFHrj34O5NXXXuWlRS9x9+i77QJ7C7FiYYxJujR/GlNvn8rC5xfW3LH12k9fqzkaUFXCp8KEqkM1xSNUHaKquopQJPTd8qj5N7e8yb7IPjgF4XCYR+c9yoKqBYzOHc2YPmPIy44xgJRpUEKLhYiMBGYDPuBlVf1trfenA1OACPAVcI+qfhr1fhawE1isqj9JZFZjTHI1dMeWiBDwBQj4AmSSGdf6etzRg+IXiwmFQgSDQab/3XRKs0p55eNXeGn7S/Tt2JcxfcYwKneU3ZnVBAkrFiLiA54HbgQ+AzaKSJGqfhLVbAtwhaqeEJFpwO+AO6LefwpYk6iMxpjW1ZJ3bNX3HMiRyiMsL1/O0rKlPLvpWZ7b9BwFOQWM7jOam3rdRIdghxbZfjyOh4/z9sq3WfynxXwR+oJhw4aR7k8nzZdGmj8Nf0rj/0tev359i98iXZdEHlkMAfaqahmAiLwBjANqioWqro5q/2dg0ukZESkAcoB3gSsSmNMYc46o6zmQzumdmdh/IhP7T2T/sf0s3beU4rJiZq6fyawNs7im+zWM6TOGEReOIN2f3iI5ToRPUPZtGXu/2UvpN6U1f5Z+VOpcVwkrb73+1ve6YUlNSSXNn0a6P72miKT702uW1X7v4I6DzL5vNpGTEYLBIKtWrUpYwUhksegOHIia/wwY2kD7e4FlACKSAjwDTAYKExXQGNO29Mzqyf0D7+e+y+5j59GdFJcVs2zfMkoOlNDO344bet3A6NzRDO02NK5v+ZWRSsq+LTuzIHxTyucVn9e0CaQEyO2Qy+VdLifjWIbTv5Y7HO+wymHcMvQWqqqrOBE5QVWkispIJVWRqprpyupKKiOVHK06+r3lh5Yc4mToJCicPHmSkpKSs7JY1HULQp23N4jIJJyjhxHuogeApap6oKE7GURkKjAVICcnh5KSkiaHraioaNbnE83r+cD7Gb2eDyxjS2hMviu5koLzC9hTtYdNJzbxXtl7FJUWkZmSyeD2g8n+NJsvPvmCSy67hE59O3H45GEOhQ9xKHyIw+HDfB35uuauLT9+uqR2oWtqVwZ1GES3QDe6pXajk78TPvHBKeiZ35OlqUsJh8Ok+lO5Kvcqcg7n1B/Q577qGYRw+4jtPLzkYSLhCH6/n6ysrMT9bFQ1IS/gamB51PxjwGN1tLsB5yJ2l6hl84H9QDlwBDgG/Lah7RUUFGhzrF69ulmfTzSv51P1fkav51O1jC2hOfmqIlW6snyl/mz1zzT/X/JVAqKkoBIQ7fNEH7103qU66NVBOv7t8frzkp/r3K1zdUX5Ci39plTD1eG4trFu3TqdMmWKrlu3rsk5a69v1qxZTV4f8KHG8X96Io8sNgIXi0gu8DlwJ3BXdAMRuRx4ARipql+eXq6qE6Pa/BDnIviMBGY1xhiCviCFvQop7FVI1gdZPBV5qqbn3mtD1/LEuCfomdWzWUPZ1nVdpTmS1a1LSqJWrKoR4CfAcpwjhz+p6g4RmSkiY91mvwcygDdFZKuIFCUqjzHGNMbNhTcTDAbx+XwEA0Gm3jaVvOy8xI557mEJfc5CVZcCS2stezJq+oY41jEPmNfS2YwxpiGJ6rn3bGVPcBtjTD2S3XOvlyXsNJQxxphzhxULY4wxMVmxMMYYE5MVC2OMMTFZsTDGGBOTFQtjjDExifO099lPRL4CPo3ZsH6dcboW8Sqv5wPvZ/R6PrCMLcHr+cBbGXup6vmxGp0zxaK5RORDVfVsV+hezwfez+j1fGAZW4LX88HZkbE2Ow1ljDEmJisWxhhjYrJi8Z0XWztADF7PB97P6PV8YBlbgtfzwdmR8Qx2zcIYY0xMdmRhjDEmpjZTLESkh4isFpGdIrJDRB6qo81EEdnmvtaJyEAPZhzn5tsqIh+KyDVeyxjV9koRqRaRCV7KJyLXici37j7cKiJP1rWu1swYlXOr22aNl/KJyCNR++9j9+d8nscydhCRJSLykdvmHz2Wr6OILHb/Pf+fiFyarHxNEs9weufCC+gGDHanM4HdwIBabYYBHd3pUcAGD2bM4LvTh5cBu7yW0X3PB7yPM57JBC/lA64D/rc1fg8bkTEb+ATo6c538VK+Wu1vBd734D58HHjanT4fOAoEPJTv98Av3el+wKpk7sPGvtrMkYWqHlLVze70X3FG7+teq806Vf2LO/tn4EIPZqxQ97cLaA8k9aJTPBldPwUWAV/W8V7CNCJfq4kz413AW6q6322XtP3YhH3498CCZGQ7Lc6MCmSKiOB8yToKRDyUbwCwym2zC+gtIjnJyNcUbaZYRBOR3sDlwIYGmt0LLEtGnro0lFFEbhORXUAxcE9yk52Rozd1ZBSR7sBtwB+Tn+qMHL2p/+d8tXt6YpmIXJLUYFEayJgPdBSREhHZJCJ3JzsbxP63IiLtgJE4XwxaRQMZ/wD0Bw4C24GHVPVUUsPRYL6PgNvdNkOAXiT5C2qjtPahTbJfON8wNgG3N9DmepxvAp28mtFtdy2w0msZgTeBq9zpeSTxNFSc+bKADHd6NLDHg/vwDzhHt+1xuobYA+R7JV9UmzuAJa2x/+LYhxOA5wABLgL2AVkeypcF/CewFXgN2AgMbK19GfPv0toBkvyDSwWWA9MbaHMZUJrsf5iNyVir/T6gs5cyupnK3VcFzqmo8V7JV0f7cg/uwxnAr6LmXwF+4JV8Ue0WA3clc981Yh8WA8Oj5t8HhnglX6224v4eJrWYNebVZk5DuectXwF2quqz9bTpCbwFTFbV3cnM524/nowXue0QkcFAAPjaSxlVNVdVe6tqb+B/gAdU9W2v5BORrlH7cAjO6VhP7UPgHWC4iPjdUz1DcY52vZIPEekAjHCzJlWcGfcDhW77HKAvUOaVfCKSLSIBd3YKsFZVjyUjX1P4WztAEv0NMBnYLiJb3WWPAz0BVPWPwJNAJ2Cu+39JRJPb2Vc8Gf8WuFtEwkAlcIe6X008lLE1xZNvAjBNRCI4+/BOr+1DVd0pIu8C24BTwMuq+rFX8rnLbgNWqOrxJOWKFk/Gp4B5IrId55v7o6qarJ5e48nXH/gvEanGufPt3iRlaxJ7gtsYY0xMbeY0lDHGmKazYmGMMSYmKxbGGGNismJhjDEmJisWxhhjYrJiYdokEflntzfQ0z34DnWX/5P7XEN9n3tZRAa40xVJyHl/a3X1YUw0u3XWtDkicjXwLHCdqoZEpDNOb6QHRaQcuKKu+/FFxKeq1VHzFaqakcCcflVNSsd3xsRiRxamLeoGHFHVEICqHnELxYPABcBqEVkNTkEQkZkisgGn88ESETnjQU0R6Swi60VkjDv/iIhsdI9afl1XAHe9z4jIZhFZJSLnu8tLRGSWOONXPCQivxKRh933LhKRlW4HiJtFJC/e7RnTXFYsTFu0AughIrtFZK6IjABQ1Tk4PZRer6rXu23bAx+r6lBV/aD2itxuJIqBJ1W1WERuAi4GhgCDgAIRubaODO2Bzao6GFgD/DLqvWxVHaGqz9T6zHzgeVUdiDP2yqFGbM+YZrFiYdocVa0ACoCpwFfAQhH5YT3Nq6m/++1UnPEIfqGq77nLbnJfW4DNOIPaXFzHZ08BC93p/waiRzxcWLuxiGQC3VV1sft3qFLVE43YnjHN0pb6hjKmhnvtoQQocfsO+gec7tRrq4q+TlFLBKf76Ztxjg7A6YPo31T1hcZGipquq68lqedzTd2eMY1iRxamzRGRviIS/e17EPCpO/1XnGEw46E4g0/1E5EZ7rLlwD0ikuFuq7uIdKnjsyk4HRqCMyre905xnbEhpzfSz0RkvLveoHvXVrzbM6ZZ7MjCtEUZwL+LSDbO0cFenFNSAC8Cy0TkUNR1i3qparWI3AksEZFjqjpXRPoD692eiyuASXx/eNnjwCUisgn4FmcQoVgmAy+IyEwgjDO+xYo4t2dMs9its8a0gkTfdmtMS7PTUMYYY2KyIwtjjDEx2ZGFMcaYmKxYGGOMicmKhTHGmJisWBhjjInJioUxxpiYrFgYY4yJ6f8BS1eZzygdGfcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#绘制牛顿法求解的波动率曲线\n",
    "plt.plot(k_16,sigma_16_newton,label='16-days',lw=1.5,)\n",
    "plt.plot(k_16,sigma_16_newton,'r.')\n",
    "plt.plot(k_51,sigma_51_newton,label='51-days',lw=1.5)\n",
    "plt.plot(k_51,sigma_51_newton,'g.')\n",
    "plt.plot(k_79,sigma_79_newton,label='79-days',lw=1.5)\n",
    "plt.plot(k_79,sigma_79_newton,'k.')\n",
    "\n",
    "plt.grid(True)\n",
    "plt.xlabel('Strike price')\n",
    "plt.ylabel('Implied volatility')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEKCAYAAADjDHn2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XlcVNX/x/HXmYFhB3dcUERxXxO3cQMkt8wsS8u1NFu1zVa/lZXtmWmlZam4oqY/Ky01MwSXQFMUd1FUVDQzURFEGGbm/P6YQdEQEBhmgPN8PO6DmTv33nkP1nw4995zjpBSoiiKoij50dg7gKIoiuL4VLFQFEVRCqSKhaIoilIgVSwURVGUAqlioSiKohRIFQtFURSlQKpYKIqiKAVSxUJRFEUpkCoWiqIoSoGc7B2gpFSrVk3Wr1+/yPtfvXoVDw+PkgtUwhw9Hzh+RkfPBypjSXD0fOBYGePi4i5IKasXuKGUslwsQUFBsjiioqKKtb+tOXo+KR0/o6Pnk1JlLAmOnk9Kx8oI7JSF+I5Vp6EURVGUAqlioSiKohRIFQtFURSlQOXmAreiKBVPdnY2ycnJZGZmXl/n4+PDoUOH7JiqYPbI6Orqip+fH87OzkXaXxULRVHKrOTkZLy8vKhfvz5CCADS0tLw8vKyc7L8lXZGKSUpKSkkJycTEBBQpGOo01CKopRZmZmZVK1a9XqhUPImhKBq1ao3tcDulCoWZUFsLPUiIiA21t5JFMXhqEJROMX9PanTUA7CZJZcuZbNlcxsUq9ZlivXjDj/tY3Q8cOon22AiAiIjAS93t5xFUWpYFSxsIGz6zZijori77s6c7pJm+tf/Km5isEVa0FIy7SsT88y5nmsZ2N/pqfBgEaakQYDIjpaFQtFcSBjxozh119/pUaNGuzfv//6+q+//poZM2bg5ORE//79+eyzz/I9TlJSEvfee+9Nx3AkqliUtNhYqg68B60xm6paJz5+5EN21WkGgKeLE96uTni7OePt5kzdKu54uzrj4+aMt5sTPm7Wx67O+LhbflYL80A7cAXGzCzMWid0ISH2/XyKotzkscceY/z48YwaNer6uqioKFatWsXevXtxcXHh/PnzdkxYMlSxKGnR0ehMRoQ0o5Um5vpfRU7shberE07aIlwiqhkCkZH8NmU24R7N+cC/Bc1LPLSiKEXVo0cPkpKSblr37bff8sYbb+Di4gJAjRo18tw3Li6OMWPG4O7uTrdu3a6vT0pKYuTIkVy9ehWAGTNm0KVLF0aOHMlDDz3EwIEDARg+fDgPP/wwDRs2ZPTo0RgMBsxmMytXrqRRo0Yl+jlVsShpISEIFx0YDAidjsr9e4OHrnjH1OsRz2RydKuBmdGJzBzWrmSyKko58t4vBzh49gomkwmtVlsix2xe25t3BrS44/2OHDnCli1bePPNN3F1deXzzz+nQ4cO/9lu9OjRfP311wQHB/Pqq69eX1+jRg02bNiAq6srR48eZejQoezcuZOxY8cybdo0Bg4cSGpqKjExMSxYsICXXnqJF154geHDh2MwGDCZTMX63Hmx6d1QQoi+QogEIUSiEOKNPF5/WgixTwgRL4TYKoRonuu11kKIWCHEAes2rrbMWmL0estF6PffL9GL0R7OglFd/Fm7728Sz6eXyDEVRbENo9HIpUuX2LZtG1OmTGHIkCFYxuy7ITU1lcuXLxMcHAzAyJEjr7+WnZ3NE088QatWrRg8eDAHDx4EIDg4mMTERM6fP8/SpUt58MEHcXJyQq/X89FHH/Hpp59y8uRJ3NzcSvwz2axlIYTQAjOBXkAysEMIsVpKeTDXZkuklLOs298HfAH0FUI4AYuBkVLKPUKIqkC2rbKWOL3eJhehx3QNIHxrEt9EJ/LFkLYlfnxFKctyWgCO0CnPz8+PQYMGIYSgY8eOaDQaLly4wGuvvcbu3bupUaMGy5cvv+3trNOmTcPX15c9e/ZgNptxdb3xt/LIkSOJiIhg2bJlhIeHAzBs2DA6derEmjVr6NOnD3PmzKFnz54l+pls2bLoCCRKKY9LKQ3AMmBg7g2klFdyPfUAckpvb2CvlHKPdbsUKWXJt6vKmKqeLgzvVI9V8Wc5lZJh7ziKotzG/fffz8aNGwHLKSmDwUC1atWYN28e8fHxrFy5kkqVKuHj48PWrVsBiIiIuL5/amoqtWrVQqPRsGjRoptOKz322GNMnz4dgBYtLAXy+PHjNGjQgOeff5777ruPvXv3lvhnsmWxqAOczvU82bruJkKIcUKIY8BnwPPW1Y0BKYRYL4TYJYR4zYY5y5QnejRAqxF8u+mYvaMoigIMHToUvV5PQkICfn5+zJ07lzFjxnD8+HFatmzJI488woIFC/JsRcybN49x48ah1+tvOnX07LPPsmDBAjp37syRI0dumijJ19eXZs2aMXr06OvrfvjhB1q2bEnbtm05fPjwTXdmlRRx63m0EjuwEIOBPlLKsdbnI4GOUsrnbrP9MOv2jwohXgHGAR2ADCASeEtKGXnLPk8CTwL4+voGLVu2rMh509PT8fT0LPL+tpY738KDWWw6bWRKsBtVXB2nE35Z+h06KpXxzvj4+BAYGHjTupK8wG0rxcmYkZFB586d2bJlCz4+Pne0b2JiIqmpqTetCw0NjZNSti9w58LMkFSUBdAD63M9nwhMzGd7DZBqffwIMD/Xa28Dr+b3fhVpprzTF6/KhhPXyHdW7bdfoDyUpd+ho1IZ78zBgwf/s+7KlSt2SHJnippxw4YNsm7dunLatGlF2j+v3xcOMFPeDqCRECJACKGzFoDVuTcQQuS+Ebg/cNT6eD3QWgjhbr3YHQzkvjBeoflVdmdQuzos/esU/6Zl2TuOoiil5O677+bUqVO8+OKLpf7eNisWUkojMB7LF/8hYLmU8oAQYrL1zieA8dZbY+OBCcCj1n0vYbkzagcQD+ySUq6xVday6JmQQLJNZuZsPW7vKIqiVAA27ZQnpVwLrL1l3aRcj1/IZ9/FWG6fVfIQUM2DAW1qszj2JE/3aEjl4nb8UxRFyYfjXB1V7ti40ECuGkzMi0mydxRFUco5VSzKsMa+XvRtUZP5f57gSmbZ6bOoKErZo4pFGTe+ZyBXMo0sij1p7yiKUmHVr1+fVq1a0bZtW9q3t9yFumLFClq0aIFGo2Hnzp2FOk5SUhItW7a0ZdQiU8WijGtZx4fQJtWZu/UEGYa858RQFMX2oqKiiI+Pv14YWrZsyY8//kiPHj3snKxkqGJRDozv2YiLVw0s2X7K3lEURbFq1qwZTZo0KXC7uLg42rRpg16vZ+bMmdfXJyUl0b17d9q1a0e7du2IiYkBLGNDrVq16vp2w4cPZ/Xq1Rw4cICOHTvStm1bWrduzdGjR//zXsWhhigvB4L8K9OlYVW+33ycEZ39cXV27N6rimIT696Ac/twMxlBW0JfbTVbQb9PCtxMCEHv3r0RQvDUU0/x5JNPFvotysow5aplUU6M7xnI+bQsVsQl2zuKolQ4f/75J7t27WLdunXMnDmTzZs3F2q/sjRMuWpZlBP6BlUJ8q/MrOhjPNKhLs5FmZVPUcoyawvgmh2GKK9duzZgaQ088MAD/PXXX7e9VjF69Gji4uLw8/NjyZIlZWaYcvWNUk4IIRjfM5Azl6/x0+4z9o6jKBXG1atXSUtLu/74999/z/eOpnnz5vHnn3+ydu3aMjVMuSoW5UhI4+q0rOPNN1GJGE1me8dRlArhn3/+oVu3brRp04aOHTvSv39/+vbty08//YSfnx+xsbH079+fPn365Ll/mRmmvDCjDZaFpSKNOpufdfv+lv6v/yp/3p1s20B5KC+/Q3tSGe9MRRt1NsfVq1dlgwYN5OXLl+9oP0cddVaxg97NfWns68mMjYmYzbaZq0RRFPv5448/aNq0Kc8999wdz2dRHOoCdzmj0QjGhQbywrJ4fj94jr4ta9k7kqIoJShnmPLSploW5dC9rWsTUM2Drzcm5kwepSiKUiyqWJRDWo3gmZCGHDh7heiEf+0dR1GUckAVi3LqgbvqUKeSG19tPKpaF4qiFJsqFuWUs1bD0yEN2X3qMrHHUuwdR1GUMk4Vi3JscJAfNbxc+Hpjor2jKEq5lZCQQNu2ba8v3t7eTJ8+nT179qDX62nVqhUDBgzgypUrBR5LDVGu2IWrs5YnezQg9ngKO5Mu2juOopRLTZo0IT4+nvj4eOLi4nB3d+eBBx5g7NixfPLJJ+zbt48HHniAKVOm2DtqsahiUc4N61SPKh46ZkSp1oWi2FpkZCQNGzbE39+fhISE6+ND9erVi5UrV+a5jxqiXHEI7jonHu8WwJT1CexLTqWVX+l14lGU0vTpX59y+OJhTCYTWm3JDNPftEpTXu/4eqG3X7ZsGUOHDgUskx+tXr2agQMHsmLFCk6fPp3nPmqIcsVhjNL74+3qxIyokv1LQ1GUGwwGA6tXr2bw4MEAhIeHM3PmTIKCgkhLS0On0/1nHzVEuZUQoi/wJaAF5kgpP7nl9aeBcYAJSAeelFIezPV6PeAg8K6U8nNbZi3PvFydeaxrAF9FHiXhXBpNapbu8M2KUhpyWgBpdhiiHGDdunW0a9cOX19fAJo2bcrvv/8OwJEjR1izZg2ghij/DyGEFpgJ9AOaA0OFEM1v2WyJlLKVlLIt8BnwxS2vTwPW2SpjRTK6S308dFpmqmsXimITS5cuvX4KCuD8+fMAmM1mPvjgA55++mlADVGel45AopTyuJTSACwDBubeQEqZ+14yD+B67zEhxP3AceCADTNWGJU9dIzQ+/Pr3rMc/zfd3nEUpVzJyMhgw4YNDBo06Pq6pUuX0rhxY5o2bUrt2rVvGk48t7IyRLktT0PVAXJf0UkGOt26kRBiHDAB0AE9res8gNeBXsArNsxYoYzt1oD5fybxbfQxpgxuY+84ilJuuLu7k5Jyc+fXF154gRdeeKHAfYOCgtizZ8/15++++y4AjRo1uql18PHHH19/nJGRcf2id46JEycyceLEon6EAtmyWOR1Iu4/405IKWcCM4UQw4C3gEeB94BpUsr0253PAxBCPAk8CZZKGx0dXeSw6enpxdrf1koqX486Gn7clUxHjxSqu5dsw7Ki/A5tSWW8Mz4+PtdnqcthMpn+s87RFCdjVFQU48aNY9y4cWg0mjs6TmZmZtH/7Qoz6UVRFkAPrM/1fCIwMZ/tNUCq9fEWIMm6XAYuAuPzez81+VHhnL2cIQP/t0bO+niRlB99JGVMTIkcV8qK8zu0JZXxzlTUyY+KqjiTH9myZbEDaCSECADOAI8Aw3JvIIRoJKXMuZ+zP3AUQErZPdc27wLpUsoZNsxaYdTycWOC50VGvf04UpoQOh1ERoJeb+9oilIkUsrb3lGk3CCLOaCozS5wSymNwHhgPXAIWC6lPCCEmCyEuM+62XghxAEhRDyW6xaP2iqPcsPQjGM4m4wIkwkMBnCQUwqKcqdcXV1JSUlRIysXQEpJSkrKTbff3imb9rOQUq4F1t6yblKuxwVe/ZFSvlvyySq2Svf0xvDZxxizs9HqdIiQEHtHUpQi8fPzIzk5mX//vTFvS2ZmZrG+FEuDPTK6urri5+dX5P3VcB8VkV7PsaWr+OXLJbQaPpB+6hSUUkY5OzsTEBBw07ro6GjuuusuOyUqnLKQ8VaqWFRQTR/ozRtnPVh3xUgfs0SjUed8FUW5PTU2VAUlhODx7g04ceEqGw+ft3ccRVEcnCoWFVi/ljWp5ePKnK3H7R1FURQHp4pFBeas1fBYl/psO36R/WdS7R1HURQHpopFBfdIx3q467SEbz1h7yiKojgwVSwqOB83Z4a0r8sve8/yz5VMe8dRFMVBqWKhMLprfYxmycLYJHtHURTFQalioeBf1YNezXyJ2H6Ka4aSnYpRUZTyQRULBYDHuwVwOSOblbuS7R1FURQHpIqFAkDHgCq0quND+J8nMJvVODuKotxMFQsFsHbS6xbA8X+vEn1EddJTFOVmqlgo193TqhY1vV2Zq26jVRTlFqpYKNfpnDSM6uLPn4kpHDx7peAdFEWpMFSxUG4yrGM93Jy1hP+pWheKotxQYLEQQuwUQowTQlQujUCKfVVy1/FQkB+r489yPk110lMUxaIwLYtHgNrADiHEMiFEH6HmMCzXRnetT7bZzOLYk/aOoiiKgyiwWEgpE6WUbwKNgSVAOHBKCPGeEKKKrQMqpa9BdU/CmtZg8fZTZGarTnqKohTymoUQojUwFZgCrAQeAq4AG20XTbGnx7s14OJVAz/tPmPvKIqiOIACZ8oTQsQBl4G5wBtSyizrS9uFEF1tGU6xn84NqtC8ljdzt57gkQ51UWceFaViK0zLYrCUMkxKuSSnUAghAgCklINsmk6xGyEEY7sHkHg+nU1H/rV3HEVR7KwwxeL/CrlOKWfubV2bGl4uqpOeoii3Pw0lhGgKtAB8hBC5WxDegKutgyn2p3PS8GiX+kxZn0DCuTSa1PSydyRFUewkv5ZFE+BeoBIwINfSDniiMAcXQvQVQiQIIRKFEG/k8frTQoh9Qoh4IcRWIURz6/peQog462txQoied/rBlJIxrGM9XJ01aiY9RangbtuykFKuAlYJIfRSytg7PbAQQgvMBHoByVj6aayWUh7MtdkSKeUs6/b3AV8AfYELwAAp5VkhREtgPVDnTjMoxVfZQ8eD7fxYEZfMq32bUM3Txd6RFEWxg9u2LIQQr1kfDhNCfHXrUohjdwQSpZTHpZQGYBkwMPcGUsrcAxB5ANK6freU8qx1/QHAVQihvqXsZEy3AAxGM4u3qU56ilJR5Xfr7CHrz51FPHYd4HSu58lAp1s3EkKMAyYAOiCv000PArtz3bKbe98ngScBfH19iY6OLmJUSE9PL9b+tmbvfG2qa5m7+SjNxRl02rxvo7V3xoI4ej5QGUuCo+eDspHxP6SUNlmAwcCcXM9HAl/ns/0wYMEt61oAx4CGBb1fUFCQLI6oqKhi7W9r9s639ei/0v/1X+UPf5267Tb2zlgQR88npcpYEhw9n5SOlRHYKQvxnZ7f3VC/YD0tdJsic18BdSgZqJvruR9w9jbbguU01be53t8P+AkYJaU8VsB7KTbWpWFVmtb0Ys7W4wxu76c66SlKBZPfaajPi3nsHUAjawe+M1gGJByWewMhRCMp5VHr0/7AUev6SsAaYKKU8s9i5lBKQM5Meq/+3162Jl6ge6Pq9o6kKEopyu9uqE3FObCU0iiEGI/lTiYtEC6lPCCEmIyl2bMaGC+EuBvIBi4Bj1p3Hw8EAm8LId62rustpVTzfdrRfW1r8+lvCczZckIVC0WpYPI7DbVcSjlECLGPPE5HSSlbF3RwKeVaYO0t6yblevzCbfb7APigoOMrpcvFScsovT9fbDjC0X/SaOSrOukpSkWRX6e8nC/ye7m5U17OolRAwzvVw8VJo2bSU5QK5rbFQkr5t/Xhs1LKk7kX4NnSiac4mqqeLgxqV4cfd50hJf0/dzMrilJOFWYgwV55rOtX0kGUsmNM1wCyjGYitp+ydxRFUUpJfj24n7Fer2gihNibazkB7C29iIqjaeTrRXDj6iyMPUmWUc2kpygVQX4tiyVYrk2s5uZrFUFSyhGlkE1xYGO7B3AhPYvV8fl1nVEUpbzI75pFqpQySUo51Hqd4hqWu6I8hRD1Si2h4pC6BVajia8Xc7eeyOltryhKOVbgNQshxAAhxFHgBLAJSALW2TiX4uByOukdPpdGzLEUe8dRFMXGCnOB+wOgM3BEShkAhAGqV7XCfW1rU81Tp2bSU5QKoDDFIltKmQJohBAaKWUU0NbGuZQywNVZy4jO/mw8fJ7E8+n2jqMoig0VplhcFkJ4ApuBCCHEl4DRtrGUsmJEZ390ThrmqU56ilKu5TeQYI6BQCbwEjAc8AEm2zKUUnZU83ThgbZ1OLZ6AyPEIXBxAb3e3rEURSlhBRYLKeXVXE8X2DCLUkY963qeGosnojMZYVkEREaqgqGUD7GxEB0NISEV/r/p/AYSTOPmAQSF9bkApJTS28bZlDLCf98OTGYjWmlGGgyI6OgK/z+WUg7ExmLqGQZZWWhcXRAV/I+g/PpZeEkpvXMtXrl/lmZIxcGFhKBxccEoNGRpnMjq1t3eiRSlWC6kZ7Hmq6XIrCy00gwGg6WFUYEV5poFQog2QM43wGYppRruQ7lBr0dERvLX13P4XNuMu1Kr8HbBeymKwzGZJUu2n2TK+gSa6PzprdMhjdkInc5yKqoCK7BYCCFeAJ4AfrSuihBCfC+l/NqmyUpRbGwsERERuLi4oK/Azcxi0esxZGXRMrUac7eeoGfTGnQNrGbvVIpSaLtPXeLtVfvZf+YKXQOr8t6zT+H8dBd1zcKqMC2Lx4FOORe6hRCfArFAuSgWMTExhIaFYjQYiYiIIDIyUhWMYpjYrxlbj17glRV7+O3FHvi4Ods7kqLk69JVA5+tP8yyHaep4eXC10Pv4t7WtSzzzNfQV/gikaMw/SwEkHtoUZN1Xbnw0/qfMBgMmM1mDAYD0RX8vGRxuem0THu4LefTsnhn1X57x1GU2zKbJcv+OkXPqdEs35nM2G4BRL4cwoA2tS2FQrlJYYrFPGC7EOJdIcS7wDZgrk1TlaJBfQeh0+lAA8JJEBwcbO9IZV6bupV4vmcjfo4/y6971ai0iuPZfyaVQd/G8MaP+2hUw4u1z3fnzf7N8XQp1GXcCqkw/Sy+EEJEA92wtChGSyl32zpYadHr9URvjOb1b17nTP0zHPI6RBe62DtWmTcutCEbE87z5k/7ae9fhZo+rvaOpCikXstm6u8JLN52kioeOr4Y0oYH7qqjWhKFUJgL3F8CP0gpvyqFPHah1+t5L/M91mnWMX3XdOp516OXf14TBCqF5aTVMG1IG+75aguv/t8eFo7pqP6HVOxGSsmPu87w8bpDXLxqYJS+Pi/1aqyuqd2BwpyG2gW8JYRIFEJMEUK0t3UoexBC8H6392lTvQ3/2/I/9l9Q59uLq0F1T97s35wtRy+waNtJe8dRKqjD564w5LtYXl6xh7pV3Fk9vhvv3tdCFYo7VGCxkFIukFLeA3QEjgCfWue3KJAQoq8QIsFaaN7I4/WnhRD7hBDxQoitQojmuV6baN0vQQjR5w4+U5G5aF34MvRLqrpV5bmNz/F3+t+l8bbl2ohO9QhuXJ2P1h5SI9MqpSotM5v3fz1I/6+2kng+nU8fbMXKp7vQso6PvaOVSYVpWeQIBJoC9YHDBW0shNACM4F+QHNgaO5iYLVEStlKStkW+Az4wrpvc+ARoAXQF/jGejybq+pWlRk9Z5BpzGTcxnFczb5a8E7KbQkhmPJQa1ydtUxYHk+2yWzvSEp5FhtL3YgItixcTdjUTYT/eYKHO9Rl48shPNyhHhqNOhVaVIWZKS+nJTEZ2I9lDu4BhTh2RyBRSnlcSmkAlmEZwfY6KeWVXE89uDEW1UBgmZQyS0p5Aki0Hq9UBFYOZGrwVI5fPs6rm17FaFYjshdHDW9XPn6gFXuTU5mxMdHecZTyKjYWc88w6s8Np/3jg+mRkshPz3blowdaUdlDZ+90ZV5hWhYnAL2Usq+Ucp6U8nIhj10HOJ3rebJ13U2EEOOEEMewtCyev5N9balLnS78r9P/2HJmC5/v/Lw037pc6teqFoPa1WFGVCK7T12ydxylHLlmMLF8x2kWfjwfs3UsJxezic+qXaRt3Ur2jlduFObW2VlFPHZe7T35nxVSzgRmCiGGAW8BjxZ2XyHEk8CTAL6+vsXqUJeenv6f/WtQg1CvUCIORWD4x0APrx5FPn5x5ZXP0RSU8e7Kkk06eHp+LJO7uOHiVLqnBMrD79AROErGM+lmok9ns/WMkWtG6F2zMY84O2M2ZiOdndjj48MVB8iZF0f5Hd4RKaVNFkAPrM/1fCIwMZ/tNUBqXtsC67G0bm77fkFBQbI4oqKi8lxvNBnl+D/Gy9YLWsstyVuK9R7Fcbt8jqQwGWMSL8j6b/wq3/xpr+0D3aK8/A7tzZ4ZM7ONclX8GTlkVoz0f/1XGfi/NfK5JbvktmMXpNlsljImRh4bO1bKmBi7ZSwMR/p3BnbKQnyn27K74g6gkRAiADiD5YL1sNwbCCEaSSlz7qzqD+Q8Xg0sEUJ8AdQGGgF/2TDrbWk1Wj7t8SmP/vYor2x6hYX9FtK4cmN7RCkX9A2rMrZbALO3nCCsmS+hTWrYO5JSBpy+mMGSv06xfMdpUq4aqFvFjdf7NmVwez+qebrc2FCv51RWFg3UeE4lLr/Jj6rkt6OU8mIBrxuFEOOxtAq0QLiU8oAQYjKWSrYaGC+EuBvIBi5hOQWFdbvlwEEs832Pk1Ka8nyjUuDu7M7XPb9m2JphjI8cz5L+S6jmpkZULaqXezdh85ELvPZ/e1n/Yg+qqIuPSh5MZsnGw+eJ2H6STUf+RQBhzXwZ0dmf7oHV1J1NpSy/lkUcN2bGq4fly1wAlYBTQEBBB5dSrgXW3rJuUq7HL+Sz74fAhwW9R2mp6VGTr8O+ZvRvo3l+4/OE9wnH1UkNYVEUrs6WwQYHztzKmz/t45vh7VTvbuW6f65k8sOO0yz76xRnUzPx9Xbh+Z6NeKRjXWr5uNk7XoWV30x5AVLKBlhaBgOklNWklFWBe7kxt0WF0qJqCz7u/jH7L+znza1vYpaqz0BRNa/tzcu9m7Bu/zl+2n3G3nEUOzObJVuPXuCZxXF0+WQjX2w4QsManswaEcTW13vyUq/GqlDYWWGuWXSQUj6d80RKuU4I8b4NMzm0sHphTAiawNS4qfjv9uf5ds8XvJOSpye6N2DjofO8s+oAHQOq4FfZ3d6RlFJ2ZeNmDi/7hQW6+qzxrE9ld2ce7xbAsI71qF/Nw97xlFwKUywuCCHeAhZjOS01AkixaSoH92iLR0m6ksTsfbPx9/ZnYODAgncqhs1bNzNv0bxyN5OfViOYOqQN/b7cwsvL97D0ic7qPHQFkXAujcjwnxn9zljamYy0cXJmyPfL6TS8L67OpTJYg3KHCtMpbyhQHfjJulS3rquwhBC82flNOtXsxLux77Lj3I4Sf49MYyaRpyIZ9e0oQsNCWTh/IT3NlXX0AAAgAElEQVTDehIbG1vi72VPdau4M2lAc7afuMjcrSfsHUexIbNZ8sfBfxg2ext9pm/m2oZIdGYjTtKMi9lI8N8HVKFwYIXplHcReEEI4SmlVCPBWTlrnJkaMpWR60byUvRLRNwTgb+3f7GOmZGdwdYzW9lwcgObkjdxzXiNtE1pSKMEM2QZsoiOji5XrQuAwUF+/HHwH6asT6B742o0relt70hKCUrLzGbFzmQWxCZxMiWDWj6uvNa3CSN6j0F773IwGECns8xzrTiswsxn0QWYA3gC9YQQbYCnpJTP2jqco/Nx8WFmz5kMWzuMcZHjiLgnAh+XOxvR8mr2VTYnb2bDyQ1sSd5CpimTKq5VuLfBvfTy74Ux0Eif1X3IzMoELbTs2NJGn8Z+hBB8PKgVfaZv5sVl8awa3xUXJ/UXZll3MuUq82OSWLEzmfQsI+3qVeLVPk3o06ImzloNEAiRkRAdbSkU5eyPoPKmMNcspgF9sHSUQ0q5Rwhhv3EvHExd77p8GfolY38fy4tRL/J9r+9x1uY/Tn6aIY3o09FsOLmBP8/8icFsoJpbNe4PvJ/e9XvTrkY7tBrrl2VtiIyMZOacmeyos4M4tzgGUJhxHMuWqp4ufPpgax5fsJNpG47yRr+m9o6kFIGUkthjKYT/eYLIw+fRCsG9rWsxumsAbfIap0mvV0WijChUD24p5elb7oO3Wwc5R9TOtx2Tu05m4paJTN42mcldJv+n30BqVipRp6PYcHIDMWdjMJqN+Lr7MqTJEHr596JtjbZoRN6XkPR6PVlZWTR3a87yhOWMaTmGet71SuOjlaqwZr4M7ViP7zYfo2fTGnQMyLdfqOJAMrNN/Lz7DPP+TCLhnzSqeOgYHxrIiM7++Hqr/kjlQWGKxWnrqSgphNBhGRn2kG1jlT33NriXk1dOMmvPLIwnjHif9qadvh1ptdP44+QfbP97O0ZppLZHbYY3HU6v+r1oVa3VbQtEXp5o9QQ/Hf2JWXtm8VH3j2z4aeznrf7NiDl2gblTImhTKxWXu8PUX54O7FxqJou2JbFk+ykuZWTTrJY3nz3Umvva1FYXq8uZwhSLp4EvsQwRngz8DoyzZaiy6tk2zxITE8OnL32KNEqEkyDgtQCa3NWEUS1G0du/N82rNi9yb+Xq7tUZ2nQo8w/M5/FWj9OwUsMS/gT25+HixHcNDfj/7xWczEbkxx8hIiNVwXAUsbHUi4gg4VIWMzJrsG7f35ikpFczX0Z3DaBzgyqqN345VZi7oS4Aw0shS5knhKDppaaW0azMIEyCgQxk6gNTS+x/oNEtR/NDwg98E/8NU0OmlsgxHU3ThF2YzUY0ZjOmLAOaqCiEKhZ2Z/4zBhkWhr/BgGH+Qi6N/IRHB/XhUX196lVVHSrLu9ueAxFCvGb9+bUQ4qtbl9KLWLbc3fNuXF1c0Wq1uOhcGHzP4BL9S6uya2VGNh/J7yd/51BKOT0bGBKCcHHBpNFi0GiZ6+SfM1S9Yid7Tl9m8acLkAbD9cmF5gZk8Pa9zVWhqCDya1nkfBPtLI0g5YVerycyMpLo6GhCQkJs0idiVItRLDm8hJnxM5kRNqPEj293ej0iMhJNVBRLdfX54IIPZ349yKR7i34KTymai1cNTFl/mGU7ThNaqwXDdDrM2QY0Ljpc7u5p73hKKbptsZBS/mL9uaD04pQPer3eph3nvHXejG4xmq92f8Wef/fQpnobm72X3ej1CL2e0VJy+teDzPszCY0QvNW/mSoYpcBkliz96xSf/55AWqaRMV0DePHu3jg91oHj4eE0GDNGXUeqYPKbz+IX8pjKNIeU8j6bJFIKZXiz4Sw+tJgZu2cwu/dse8exGSEEk+5tjpQwd+sJNAL+d48qGLa069QlJq3az/4zV+gUUIXJA1vSpKaX5UU1uVCFld9pqM9LLYVyx9yd3Xm85eNM2TmFHed20KFmB3tHshkhBO8MaI5ZSmZvOYFGCN7o11QVjBJ2IT2LT9cdZkVcMr7eLnw19C4GtK6lfs8KkP9pqE05j639K5piaWkkSCkNpZBNKcCQJkNYcGABM3bPYH7f+eX6f2ohBO/d1wKzlHy3+ThCCF7v26Rcf+bSYjSZidh+iqm/J5BhMPFUjwY8F9YITxdbzrqslDWFGRuqPzALOIZlprwAIcRTUsp1tg6n5M/VyZUnWz/JB9s/IOZsDF3rdLV3JJsSQjD5vpZICbM2HUMj4NU+qmAUx46ki0xadYBDf1+hW2A13r2vBYE1PO0dS3FAhfnTYSoQKqVMBBBCNATWAKpYOIBBjQYRvj+cr3d/TZfaXcr9F6dGI3h/YEvMEr6JPoZGCF7u3bjcf+6Sdj4tk0/WHubH3Weo7ePKN8Pb0a9lTfV7VG6rMMXifE6hsDoOnLdRHuUOOWudebrN00yKmcTG0xsJqxdm70g2p9EIPry/JVJKZkQlotEIJvRqbO9YZUK2yczC2JNM33CETKOJZ0MaMr5nIO46dcpJyV9h/gs5IIRYCyzHcs1iMLBDCDEIQEpZIefjdiQDGg5g7v65zNg9g9C6oXc03lRZpdEIPnqgFWYp+SryKBoBL96tCkZ+th1P4Z1VB0j4J43gxtV5Z0BzGlRXp5yUwilMsXAF/gGCrc//BaoAA7AUD1Us7MxJ48SzbZ7l9S2vsz5pPf0C+tk7UqnQaASfDGqNWcL0P44iELxwdyN7x3I4/1zJ5MM1h1i95yx1Krnx3cggejf3VaeclDtSmLGhRhf14EKIvlgGIdQCc6SUn9zy+gRgLJbRlP4FxkgpT1pf+wzoj2VIkg3AC1KN+XBbfQP6MnvfbL6J/4Ze/r1w0lSM0woajeDTB1sjJUz74wgaAc+FqYIhpeRy5GaOLv+F6cY67KzdlOfDGvFMcEPcdGo0WOXOFeZuqADgOaB+7u0L6pQnhNACM4FeWEar3SGEWC2lPJhrs91AeyllhhDiGeAz4GHrkOhdgdbW7bZiadlEF+5jVTwaoWH8XeN5MepFfj3+K/cH3m/vSKVGqxF89lBrpJRM3XAEjUYwLjTQ3rFKRXqWkaQLVzl+4SrH/03nxIWrnLhwFe9dO5i96A3amYzMc3bm4qq11FLXdZRiKMyfnz8Dc4FfAPMdHLsjkCilPA4ghFgGDASuFwspZVSu7bcBI3JewnL6S4fldl1nLKfClHz0rNuT5lWbM2vPLPoH9C9wxr7yRKsRTBncBrOUTFmfgBDwbEj5KBjZJjOnLmZw4l9LIchdGM6nZV3fTgio7eNGg+oePJp9EhezEY0042QyUmv3duirxnJSiq4wxSJTSlmUUWbrAKdzPU8GOuWz/eNYb8eVUsYKIaKAv7EUixlSynI6xGrJEULw3F3P8cwfz/BT4k8MaTLE3pFKlVYjmDqkLRL47LcENELwdHDZmfMjY9MWXGcvYn3iRf6q2fR6K+HUxQxM5htnYKt46Aio5kGPxtVpUN2DBtU8CKjmiX9V9xsTDjU3w8rvwWAAnc4yx7WiFIMo6DKAEGIY0AjLpEfX/4yRUu4qYL/BQB8p5Vjr85FARynlc3lsOwIYDwRLKbOEEIFYrnU8bN1kA/C6lHLzLfs9CTwJ4OvrG7Rs2bJ8P0t+0tPT8fR03DtDCptPSsn0f6aTYkxhUu1J6DS6Ukhn4Si/Q5NZ8v3eLLafM/FwEx39AiwtLEfJlxfDzn30mPgKTiYj2VonHh36IX83aUFND0FNdw01PQS+Hhpqumvw1BXuwrT3gQNUio/nctu2XGnRosSyOvLvERw/HzhWxtDQ0DgpZfsCN5RS5rsAH2NpFWwCoqzLxkLspwfW53o+EZiYx3Z3YxkOvUauda8Cb+d6Pgl4Lb/3CwoKksURFRVVrP1t7U7y/fX3X7Ll/JZywf4FtguUB0f6HWYbTfLZiDjp//qvcvbmY1JKx8qXW9zJi/LLsNEyW2ikBGnWaqXpww/tHeu2HPX3mMPR80npWBmBnbKA73Mp5e0nP8rlAaCBlDJYShlqXQpz8nMH0EgIEWAdW+oRYHXuDYQQdwHfAfdJKXN39DsFBAshnIQQzlgubqvTUIXUoWYHOtXqxNz9c8nIzrB3HLtw0mr48uG29G9Viw/WHGLu1hP2jpSnqMPnGTZ7G4ebtkPj4oJZo0HodGhCQ+0dTVFuUphisQeodKcHllIasZxaWo/li365lPKAEGKyECLnTqopgCewQggRL4TIKSb/h2Usqn3W998jrfNrKIXz3F3PcTHzIksOL7F3FLtx0mqY/khb+rWsyZrvVpI1YyHExto71nUr45IZu3AngTU8ee/jJ9BsjCRpzBhQc44rDqgwF7h9gcNCiB3cfM2iwPkspJRrgbW3rJuU6/Hdt9nPBDxViGzKbbSp3oZgv2DC94czpMkQvHXe9o5kF85aDV/Xz8S8/C002dkYflmCecMGXHt0t1smKSXfbz7Ox+sO0zWwKt+NbG8Z4VXNFaE4sMK0LN7BcirqIyyDCuYsioMb13YcaYY0Fh1cZO8oduW0ZTPOJiNO0owmO5tFnyxkb/Jlu2QxmyUfrDnEx+sOc2/rWoQ/1kENBa6UCQUWCynlpryW0ginFE+zqs3o5d+LRQcXcSnzkr3j2E9ICEKns1wPcNGxvV5rBn0Tw8yoxJtuSbU1g9HMS8vjmbv1BI91qc9Xj9yFi5PqTa2UDbctFkKINCHElTyWNCHEldIMqRTduLbjyMjOYN6BefaOYj96PURargdoN27k8y+eok+LmkxZn8DQ2ds4c/mazSNczTLy+IIdrIo/y6t9mvDOgOZoNGpsJqXsuG2xkFJ6SSm981i8pJQV8wR4GdSwUkP6N+jP0kNLuXDtgr3j2I9ez6nhw0Gvp5K7jhnD7uLzwW04cCaVvtM388ueszZ765T0LIbN3kbMsRQ+e7A140ID1SB+SplT/seyVnimzTNkm7OZs2+OvaM4DCEEDwX5sfaF7jSs7slzS3czYXk8aZnZJfo+py9m8NCsWA6fS+O7EUEM6VC3RI+vKKVFFYsKoJ53Pe4PvJ/lCcv5O/1ve8dxKP5VPVjxtJ7nwxrx8+4z3PPVFuJOXiyRYx/6+woPfhvDxasGIsZ24u7mviVyXEWxB1UsKoinWlvuRP5u73d2TuJ4nLUaJvRqzPKn9EgJg2fFMm3DEYymOxk382bbjqcwZFYsWo1gxdN62tevUoKJFaX0qWJRQdTyrMXgxoP5OfFnTl05VeLHj42NJSIiglgH6vR2p9rXr8LaF7pzf9s6fBl5lCHfxXIq5c57wP+2/29Ghf+Fr48rK5/pQmNfLxukVZTSpYpFBTK21VicNc7M2jOr2MeSUpJyLYW9/+5l6v9NJbhnMOHh4YSFhZXpguHt6swXD7fly0facvR8Ov2+3Mz/xSXnjFFWoMXbTvJsxC5a1PZmxVN6aldys3FiRSkdqjdQBVLdvTpDmw5l/oH5PN7qcRpWyn/47ozsDM6kn7m+JKclk5yeTHJaMmfSz3DNaLnl9N9f/yU7KxskZBmyiI6ORl/GeyEPbFuHIP/KTFi+h1dW7CEq4Twf3d8KH/e85wiRUvJl5FGm/3GUnk1rMHNYOzUjnVKuqGJRwYxuOZofEn7grYi3aJnakladWuHb3PemYpDz+GLmzRd63Z3cqeNVh7pedelcqzN+Xn74efpxvuZ5Rq0ZRWZWJlIruVj7IlLKMn97qF9ld5Y+0ZlZm44xbcMRdp28xBdD2qJvWPWm7UxmyaRV+4nYfoqHgvz4eFArnLWq0a6UL6pYVDCVXSvTLbsb0yZMY5lxGcJJEPBaAO6B7miFlpoeNfHz8iO0bih+Xn7U8ayDn6cfdbzqUNmlct4FoC74RfoxO3w2piATv/EbYovg/a7v46J1Kf0PWYK01ilauwVW48Uf4hk2ZxtP9WjIhF6N0TlpyMw28eKyeH47cI5nQhryWp8mZb5IKkpeVLGogLySvcAImEGYBH1MfXjzwTfxdffFSVO0/yT0ej1ZWVkEBwczd/9cvtz1JWfTzzI9dDrV3KqV7AewgzZ1K/Hrc914/9eDzNp0jK2J//JF3WvELfiJ854NefupB3m8W4C9YyqKzai2cgXUJ6wPri6uaLVaXHQujBwwkjqedYpcKHITQjC21Vi+CPmChIsJDF8znKOXjpZAavvzcHHikwdbM2tEEFX2xlF38AAGr/qe5Sve5nHtOXvHUxSbUsWiAtLr9URGRvL+++8TGRlpk4vRvfx7Mb/ffLLN2YxcN5LNyZsL3imX2NOxfLzlY2JPO96dVX1b1mRm7TR0ZstItk7GbIiOtncsRbEpVSwqKL1ez8SJE21611KLqi1Y0n8J9bzq8dzG54g4FFGoW1BjT8cStiCEtze+SdjCng5ZMLz63o3WxQW0WtDpICTE3pEUxaZUsVBsqqZHTeb3nU+IXwif/PUJH27/kGxz/uMvRSdFYzBlY0JiyM4kesVw2Pw5nD8MhezvYHPWkWx5/301s51SIagL3IrNuTu7My10GtN3TWfe/nmcunKKz0M+v+3sfSH1Q9A5uWIwGdBpNIS4VIaN71uWKg2haX9oei/4dQCNHf/e0etVkVAqDFUslFKhERomBE0gwDuAydsmM2LtCGb2nEld7/+OwqqvqydyVCTRSdGE1A9BX1cPV85Cwlo4vAa2fQMxX4FHDWjSz1I4AnqAs6sdPpmiVAyqWCil6oFGD+Dn5cdL0S8xbO0wpodOJ8g36D/b6evqLUUih3dt6DDWsly7DIl/wOFfYf9K2LUAdJ4QeLelcDTqBW6VSvFTKUr5p65ZKKWuQ80ORNwTQSWXSoz9fSyrj62+swO4VYJWD8Hg+fDacRj+f5bnJ2Pgx7EwpSEsvB/+mg2pZ2zyGRSlolEtC8Uu/L39WXzPYl6Ofpk3t77JidQTPHfXc2jEHf794uRiaUk06gX9p8GZnZYWx+E1sPYVy1K7HbG1WrD8n39waaBDX6+LbT6UopRjqlgoduPj4sO3vb7lo+0fMWffHJJSk/io+0e4ORVxpFaNBup2tCy9JsO/R+Dwr8TuXUJY3EwMwPx5K4lsOBB966EQEAxeakIiRSkMm56GEkL0FUIkCCEShRBv5PH6BCHEQSHEXiFEpBDCP9dr9YQQvwshDlm3qW/LrIp9OGucmdR5Eq+2f5XIU5E89ttjnM84XzIHr94Yuk9gQ8v7MJ3TQYyWrNOS6JPR8OMTMLUxfNMFfvsfHN0Ahqsl876KUg7ZrGUhhNACM4FeQDKwQwixWkp5MNdmu4H2UsoMIcQzwGfAw9bXFgIfSik3CCE8gaJPW6Y4NCEEo1qMwt/bn9c2v8bQNUMZ6zWWpN1JhISE/KfjoJSS9Ox0Uq6lkJKZku/Pi5kXSTmcQva8bGS2RDgJXCPehq6hcDwajkXBjjmwbSZonKFuJ2gYAg1CofZdoFHDjCsK2PY0VEcgUUp5HEAIsQwYCFwvFlLKqFzbbwNGWLdtDjhJKTdYt0u3YU7FQQTXDWZhv4WM+HYEwycPByM4OTsx4qsRuDd0v6kIGMyG/+wvEFR2rUwV1ypUdatKm+ptqOpWle1/beeE8QRIkCbJlyu/olFQY/p3fRHR7SXIvganYi2F43g0bPzAsrj6WG7JbRBiKR5VGoAaUVapoERhZwC74wML8RDQV0o51vp8JNBJSjn+NtvPAM5JKT8QQtwPjAUMQADwB/CGlNJ0yz5PAk8C+Pr6Bi1btqzIedPT0/H09Czy/rbm6Pmg5DKGLwpn0fxFlrakBuo9WI9m9zfDS+uFl9YLb603Xhqv689z1nloPNCK/7YEDhw4wMsvv0x2djZOzk50mNiB1PqpNHdrziNVHqGyU+Wbtnc2pFL50p7ri2vWvwBkutTgYpU2XKrchsuV2hB/7TTxqfG09WlLC58Wxf7cULH+nW3F0fOBY2UMDQ2Nk1K2L3BDKaVNFmAwMCfX85HA17fZdgSWloWL9flDQCrQAEvrZyXweH7vFxQUJIsjKiqqWPvbmqPnk7LkMsbExEg3Nzep1Wqlm5ubjImJKZFjjh07VsbExEijySgXH1wsOyzuIDsu7iiXHloqTWZT3juazVJeSJRy+/dSLh0m5Ud1pXzHW8a84y7d3tVI7btCur3vImOSNhc7o5QV69/ZVhw9n5SOlRHYKQvxnW7L01DJQO7uuX7A2Vs3EkLcDbwJBEsps3Ltu1veOIX1M9AZmGvDvIqDyBkVNzo6Os9rFkU9ZlZW1vVjDW82nJC6IbwX8x4fbv+QdSfW8W6XdwnwuWVOCiGgakPL0vEJMBnh73iiN76N4cRvmACDMYvoRQPQNx4IjXpbOgequ6yUcsaWxWIH0EgIEQCcAR4BhuXeQAhxF/AdltNV52/Zt7IQorqU8l+gJ7DThlkVB6PX620+j3cdzzp81+s7Vh9bzWc7PuOh1Q/xTNtneLTFozhr8p5rG60T+LUnJHQSutObLONXabWENAiDU9vg4M+W7Wq1sRaOXuDXXl0oV8o8mxULKaVRCDEeWA9ogXAp5QEhxGQszZ7VwBTAE1hhnYrylJTyPimlSQjxChApLC/EAbNtlVWpuIQQDAwcSNc6Xflo+0d8uetL1iet570u79G8avPb7pfn+FVSwrl9cPR3y3AkW6bC5ingVhka9rQUj4Zh4Fm9FD+hopQMm3bKk1KuBdbesm5Srsd357PvBqC17dIpyg3V3KrxRcgXRJ6M5IPtHzBszTAebfEoz7R5BlenvAco/M/4VUJArdaWpccrcO2S5Q6roxssxWP/SkBYbslt1MtSPNTtuUoZoXpwK0ouYf5htK/Zni/iviB8fziRpyJ5V/8u7WsWfLPIf7hVhpaDLIvZDOf2wNE/LC2PzVNg06fgVsVyjaNRL2LdvIk4tRKX0y43FyFFcQCqWCjKLXxcfHivy3v0C+jHuzHvMnr9aIY0HsJLQS/hqSvi7Y4ajaUVUfsuCH4VMi7CsY3XWx2x+5YQRgYGARFJ4UT690Xv2xo8fcGr5s0/3Sqr/h5KqVPFQlFuo3Otzvx434/MjJ/J4kOL2ZS8iUn6SfTw61H8g7tXsYyU2+ohMJuJXvc8hp3fWGYHlGaiz+1GfyYeDHn0R9W6WIuHb64iUtP6/MbP2EtHiT655cY1FUUpBlUsFCUf7s7uvNrhVfrU78M7Me8wLnIc/QL68UbHN6jiWqVk3kSjIaT1cHTx4WQZs9A5uRAybCXU1UNWOqT/A2nnIP0cpP1z88+UY3DyT8v1kVxiMRImMjAAOqElMmg8+hYPQa224OIYncHAMt/6TTcJlMDxIk5FlMypPCmJPRFNdNJGQhr2Ru/fvdj5yjJVLBSlEFpXb83ye5czZ/8cvt/7PbFnYxnoNJD0w+mEBIfQvWt3tMW4UK2vq2d6i+nMjpjNE8OfuPFF5+JpWao2zP8AxixrUbEUkej4eRiO/GxtqZiI3jEL/Y55IDRQvRnUaWe5pbdOkOW5tpS+CrLSICURLhwl9ngkYXvnYpAmS0Gr3w+9l58li1ZnWTQ5j52ti84yhlfO45z1GmdiLx0jLGoiWSYDEacWEtn1bfTefpbhXAxXLT+zr0F2Rq7lGhhyPc6+BtmWbWOzr1hODQK6mClEPhpVoVtoqlgoSiE5a515ps0z9KrXi2fnPstrb76GzJZMdp5MwGsBeDbyRKfR4axxxlnrjJPGyfJY44xOq7v+2FlrXafRXd/un4P/sOTFJZiyTexfs59Wka3urJ+JkwtUqmdZgBDPauiO/2btB6IjZMiPgBaSd8KZOMt8H7sXWT+Yu6XFkbuA+NQt+nURswlST8OFREg5CheOwIWjliKR9vf1zaKFAQMmS8dGaSL67zj0KUlgzgaTAUzZ1sUAN4/0k6dosjCQhVmAwWggetMH6HG5sYHQgLMHOLtZF3fQuVt+etawrrvxevS5HRiS/rAUXLOR6KRoVSwURSm8wMqBBBuCiTRGIqVEmAStUlsR1iqMbHM22aZsy89cjw0mw4115mzSs9Nv2i5hSwLGbCOYITMrkwlzJvB+7fcJ9gvGWXubDoL5yLMfCFhu2QVLn5BLJyA5zlI8zsRZZhaMnWF53aOGpWjUCQK/IKjdDtwq3Xyap3rzXAXhaK6fx8CUdSOMqw9Ua2wZjLFaIFRtBNUaEZJxHl1EvxsFLefUW17M5ryLiPnG45CzcejWPo3BlI3OyZmQ+xdBva43CoNWd0cFMOR0LLqFW2/kqx9yx/8O5YkqFopSBD1De/Khy4cYDAZ0Oh2vD30d/V1F/6sztmYsYT+FkZWVhdZZiyHAwEvRL1HJpRL3BNzD/YH307RKU8QdfNn9px9IbkJYRtGt0gBaD7asMxrg/AFr62OXZdbBI+tuZPSuRVh6IgZpIuLEXCKlG/qcrxChhcr1oVojCAy7XhCo2gg8quX5Ja2nWd4FLS8aDWhcLC2o233emq2IrN6E8KhwxoSOKXYr4LYFt4JSxUJRiqCkx6/KOV54eDhjxoyhQ6cObPt7Gz8n/syKIytYcngJjSs3ZmDDgfRv0J+qblVL6JPk4qS7cXtvjsxUa+GII3rfIgxp1tNGQHRgCPqgZyythsr1LfvfoXwLWhHo6+rJqpdVYscs6XxlmSoWilJEJT1+1a2DHXar041udbqRmpXKbyd+Y9WxVUzZOYVpcdPo7tedgYED6eHX4/bjWJUEVx9oGAoNQwkJ6IpuYdiNO7aC37r9aSOl3FHFQlEcnI+LDw83fZiHmz5M4qVEVh9bzS/HfyHqdBRVXKtcP03VpEoTm+bIOS1TUqd5lLJFFQtFKUMCKwcyof0Enm/3PDFnY/g58Wd+SPiBxYcW07RKU+4PvJ97Au6hsmvlgg9WSGZpJiM7gzRDGmcPn0VukdCImycgUMo9VSwUpQxy0jjRw68HPfx6cDnzMmtPrGXVsVV88tcnfL7zc4L9gmmU2ohLBy/RvUd3WgS1IM2QRrohnTRDGmnZNz++9f3FTMkAAA6fSURBVLXcz9Oz05FIMhIzOPHZCaRRsnDRQpb+spQHez1o71+FUkpUsVCUMq6SayWGNRvGsGbDOHLpCKsSV7HktyV89cFXSKPk/9u79+CoyjSP498n3elOIAlBkIBoIImGi5YgERBWRCteuCigC6sr4M4qAzjOLKujI7quM4M17DhTauGO7HirxXVZvCzikg1MECRQFgyLXAQRR0yIqICKjDJA0nSHZ/84h9jEJN25dOdAnk9VF+ecfvucHyeBp8/tfcUv5P0sj04Xdmrw8ymSQufUzmQFsshIzSAzkMl5GefVzWcEMsgKZLFqxyqqIlXoSSUcDnPPc/ew9MRSJhRMYEzfMWSnZSf5b26SyYqFMWeRwq6FPDD0AUJvhdhRuwM96TwHMvTYUKaNmnZaQcgMZJKRmkGn1E6kSErsdU8uZNm/LSMUChEIBvjhpB+yq3YXv9r0Kx7f/Dijzx/NhIIJjOo9qkXPhhhvs2JhzFmo+Jpi5gfm1z0HMmfKHEbkt+3tvSNGjEBV+dOf/8TyiuWUVpayZt8asoPZjM0by8SCiQzsNrBZz4YY77JiYcxZKBHjmJ9ab/TtvSJC/3P60/+c/txbdC8b929kecVyln60lCUfLiG/Sz43FdzEjfk30rNzzzbJYNqHFQtjzlLJGMc8WmpKat1F9yMnjlBWVUZJRQkLti7g6a1PM7zXcCYUTKA4t5hOqQ1fPzHeZcXCGNPmsgJZTCmcwpTCKew7so+SyhJKKkp4+J2HSfenc12f65hYMJHw3jDr161v06MfkxhWLIwxCZWblcs9g+/h7kF3s+3LbZRUlFBWVcYrZa9Q9ZsqNKL4Un1MXzCdgkEFBH1Bgv4gab40gr4gaf40Ar5A3XzQH3SWR82n+Zw22zdvZ/HixQSDQSs+bcyKhTEmKVIkhaKcIopyipg7bC6zH5rN3she9KQSCUdYv249e7rsIVQboiZSg6LNWn/0cyAvvfwSzy99njvG3WEX2NuIFQtjTNKl+dOYdcssXnvmtbo7tl7+yct1RwOqSvhkmFBtqK54hGpD1NTWEIqEvlseNf/6ttfZG9kLJyEcDvPgogdZUrOEcXnjGJ8/noLsGANImSYltFiIyBhgAeADXlDVX9d7/z5gBhABvgLuVNVPot7PAnYDy1T1x4nMaoxJrqbu2BIRAr4AAV+ATDLjWt8Ft15A6XOlhEIhgsEg9/3NfVRkVfDi+y/y/M7n6de1H+PzxzM2b6zdmdUCCSsWIuIDngGuAz4DNovIclX9IKrZNuByVT0uIncDvwFujXr/MWBdojIaY9pXW96x1dBzIACHqg9RVlXGisoVPLnlSZ7a8hRFOUWMyx/H9X2up0uwS5tsPx7Hwsd4c/WbLHttGV+EvmDkyJGk+9NJ86WR5k/Dn9L8/5I3btzY5rdINySRRxbDgI9VtRJARF4BJgJ1xUJV10a1/yMw7dSMiBQBOcAfgMsTmNMYc5ao/xwIQPf07kwdMJWpA6ay78g+VuxdQWllKfM2zmP+pvlc2ftKxuePZ/T5o0n3p7dJjuPh41R+W8nH33xMxTcVdX9WvFfhXFcJK2/81xvf64YlNSWVNH8a6f70uiKS7k+vW1b/vf279rNg1gIiJyIEg0HWrFmTsIKRyGLRG/g0av4zYHgT7e8CVgKISArwBDAdKE5UQGNMx5KblcvsQbOZdeksdh/eTWllKSv3rqT803I6+TtxbZ9rGZc3juG9hsf1Lb86Uk3lt5WnF4RvKvj86Od1bQIpAfK65HFZj8vIOJLh9K/lDsc7snokNw6/kZraGo5HjlMTqaE6Uk1NpKZuurq2mupINYdrDn9v+YGSA5wInQCFEydOUF5efkYWi4ZuQWjw9gYRmYZz9DDaXfQjYIWqftrUnQwiMhOYCZCTk0N5eXmLwx49erRVn080r+cD72f0ej6wjG2hOfmGMpSic4vYU7OHLce38FblWyyvWE5mSiZDOg8h+5NsvvjgCy6+9GK69evGwRMHORA+wIHwAQ6GD/J15Ou6u7b8+OmR2oOeqT0Z3GUwvQK96JXai27+bvjEBychtzCXFakrCIfDpPpTuSLvCnIO5jQe0Oe+GhmEcOfondxfcj+RcAS/309WVlbifjaqmpAXMAIoi5p/CHiogXbX4lzE7hG1bDGwD6gCDgFHgF83tb2ioiJtjbVr17bq84nm9Xyq3s/o9XyqlrEttCZfTaRGV1et1nvX3quF/1yoEhAlBZWAaP4j+XrJokt08EuDddKbk/Sn5T/VhdsX6qqqVVrxTYWGa8NxbWPDhg06Y8YM3bBhQ4tz1l/f/PnzW7w+4F2N4//0RB5ZbAYuEpE84HPgNuD26AYichnwLDBGVb88tVxVp0a1+QHORfC5CcxqjDEEfUGK+xRT3KeYrHeyeCzyWF3PvVeFruKRiY+Qm5XbqqFsG7qu0hrJ6tYldr/ELaSqEeDHQBnOkcNrqrpLROaJyAS32W+BDOB1EdkuIssTlccYY5rjhuIbCAaD+Hw+goEgM2+eSUF2QWLHPPewhD5noaorgBX1lj0aNX1tHOtYBCxq62zGGNOURPXce6ayJ7iNMaYRye6518sSdhrKGGPM2cOKhTHGmJisWBhjjInJioUxxpiYrFgYY4yJyYqFMcaYmMR52vvMJyJfAZ/EbNi47jhdi3iV1/OB9zN6PR9Yxrbg9XzgrYx9VPXcWI3OmmLRWiLyrqp6tit0r+cD72f0ej6wjG3B6/ngzMhYn52GMsYYE5MVC2OMMTFZsfjOc+0dIAav5wPvZ/R6PrCMbcHr+eDMyHgau2ZhjDEmJjuyMMYYE1OHKRYicoGIrBWR3SKyS0TmNNBmqojscF8bRGSQBzNOdPNtF5F3ReRKr2WMajtURGpFZLKX8onI1SLyrbsPt4vIow2tqz0zRuXc7rZZ56V8IvJA1P573/05n+OxjF1EpERE3nPb/L3H8nUVkWXuv+f/E5FLkpWvReIZTu9seAG9gCHudCbwETCwXpuRQFd3eiywyYMZM/ju9OGlwIdey+i+5wPexhnPZLKX8gFXA//bHr+HzciYDXwA5LrzPbyUr177m4C3PbgPHwYed6fPBQ4DAQ/l+y3wc3e6P7Ammfuwua8Oc2ShqgdUdas7/Rec0ft612uzQVX/7M7+ETjfgxmPqvvbBXQGknrRKZ6Mrp8AS4EvG3gvYZqRr93EmfF24A1V3ee2S9p+bME+/FtgSTKynRJnRgUyRURwvmQdBiIeyjcQWOO2+RDoKyI5ycjXEh2mWEQTkb7AZcCmJprdBaxMRp6GNJVRRG4WkQ+BUuDO5CY7LUdfGsgoIr2Bm4HfJz/VaTn60vjPeYR7emKliFyc1GBRmshYCHQVkXIR2SIidyQ7G8T+tyIinYAxOF8M2kUTGX8HDAD2AzuBOap6MqnhaDLfe8AtbpthQB+S/AW1Wdr70CbZL5xvGFuAW5pocw3ON4FuXs3otrsKWO21jMDrwBXu9CKSeBoqznxZQIY7PQ7Y48F9+Duco9vOOF1D7AEKvZIvqs2tQEl77L849uFk4ClAgAuBvUCWh/JlAf8ObAdeBjYDg9prX8b8u7R3gCT/4FKBMuC+JtpcClQk+x9mczLWa78X6O6ljG6mKvd1FOdU1CSv5GugfZUH9+Fc4BdR8y8CU7ySL6rdMuD2ZO67ZuzDUmBU1PzbwDCv5KvXVtzfw6QWs+a8OsxpKPe85YvAblV9spE2ucAbwHRV/SiZ+dztx5PxQrcdIjIECABfeymjquapal9V7Qv8N/AjVX3TK/lEpGfUPhyGczrWU/sQ+B9glIj43VM9w3GOdr2SDxHpAox2syZVnBn3AcVu+xygH1DplXwiki0iAXd2BrBeVY8kI19L+Ns7QBL9FTAd2Cki291lDwO5AKr6e+BRoBuw0P2/JKLJ7ewrnox/DdwhImGgGrhV3a8mHsrYnuLJNxm4W0QiOPvwNq/tQ1XdLSJ/AHYAJ4EXVPV9r+Rzl90MrFLVY0nKFS2ejI8Bi0RkJ8439wdVNVk9vcaTbwDwHyJSi3Pn211JytYi9gS3McaYmDrMaShjjDEtZ8XCGGNMTFYsjDHGxGTFwhhjTExWLIwxxsRkxcJ0SCLyT25voKd68B3uLv9H97mGxj73gogMdKePJiHn7Pbq6sOYaHbrrOlwRGQE8CRwtaqGRKQ7Tm+k+0WkCri8ofvxRcSnqrVR80dVNSOBOf2qmpSO74yJxY4sTEfUCzikqiEAVT3kFop/AM4D1orIWnAKgojME5FNOJ0PlovIaQ9qikh3EdkoIuPd+QdEZLN71PLLhgK4631CRLaKyBoROdddXi4i88UZv2KOiPxCRO5337tQRFa7HSBuFZGCeLdnTGtZsTAd0SrgAhH5SEQWishoAFV9GqeH0mtU9Rq3bWfgfVUdrqrv1F+R241EKfCoqpaKyPXARcAwYDBQJCJXNZChM7BVVYcA64CfR72XraqjVfWJep9ZDDyjqoNwxl450IztGdMqVixMh6OqR4EiYCbwFfCqiPygkea1NN79dirOeAQ/U9W33GXXu69twFacQW0uauCzJ4FX3en/BKJHPHy1fmMRyQR6q+oy9+9Qo6rHm7E9Y1qlI/UNZUwd99pDOVDu9h30dzjdqddXE32dop4ITvfTN+AcHYDTB9G/qOqzzY0UNd1QX0vSyOdauj1jmsWOLEyHIyL9RCT62/dg4BN3+i84w2DGQ3EGn+ovInPdZWXAnSKS4W6rt4j0aOCzKTgdGoIzKt73TnGdtiGnN9LPRGSSu96ge9dWvNszplXsyMJ0RBnAv4pINs7Rwcc4p6QAngNWisiBqOsWjVLVWhG5DSgRkSOqulBEBgAb3Z6LjwLT+P7wsseAi0VkC/AtziBCsUwHnhWReUAYZ3yLVXFuz5hWsVtnjWkHib7t1pi2ZqehjDHGxGRHFsYYY2KyIwtjjDExWbEwxhgTkxULY4wxMVmxMMYYE5MVC2OMMTFZsTDGGBPT/wP2dJyh7KjCegAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#绘制二分法求解的波动率曲线\n",
    "plt.plot(k_16,sigma_16_dichotomy,label='16-days',lw=1.5,)\n",
    "plt.plot(k_16,sigma_16_dichotomy,'r.')\n",
    "plt.plot(k_51,sigma_51_dichotomy,label='51-days',lw=1.5)\n",
    "plt.plot(k_51,sigma_51_dichotomy,'g.')\n",
    "plt.plot(k_79,sigma_79_dichotomy,label='79-days',lw=1.5)\n",
    "plt.plot(k_79,sigma_79_dichotomy,'k.')\n",
    "\n",
    "plt.grid(True)\n",
    "plt.xlabel('Strike price')\n",
    "plt.ylabel('Implied volatility')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.2848400904069232, 0.2842367078839068, 0.2827011409660353, 0.2790832314148648, 0.2797722610218441, 0.2751046357764691]\n",
      "[0.28484009113162756, 0.28423671051859856, 0.2827011412009597, 0.2790832305327058, 0.2797722602263093, 0.27510463539510965]\n",
      "###############\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x24c8372d710>]"
      ]
     },
     "execution_count": 110,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VNX9//HXJ5MNyEKAsAYSlrAECNu4BBUBxaJ+RRCt2qq0xaW29tvW2m+1arUuda1arVq3Vrtp3VCsUtQIiIotCVvY9yXsyBYI2c/vjxn4xTSQAZLcmeT9fDzmkdx7z535nEcg75xz75wx5xwiIiJRXhcgIiLhQYEgIiKAAkFERIIUCCIiAigQREQkSIEgIiKAAkFERIIUCCIiAigQREQkKNrrAo5Hu3btXEZGhtdliIhElPz8/F3OudS62kVUIGRkZJCXl+d1GSIiEcXMNoTSTlNGIiICKBBERCRIgSAiIoACQUREghQIIiICKBBERCSoeQTCnDnwwAOBryIiUquIeh/CCZkzh4pRo4kqL8PFxlE+fTrxI87yuioRkbDT9ANh5kysrIwoV0VFaSlP/epFPrq4igFdksnukszAtGSyOiXTItbndaUiIp5q+oEwciS++DhcWRkWE0PahAtIS27Jpyt38fa8zQBEGWS2T2RgWjIDj4REEvExCgkRaT7MOed1DSHz+/3uhJaumDMHZs6EkSMhJwcA5xzb9pdQULiPgs3BR+E+vjpYBoAvyshsn8DALslkpyUzMK01fTsmKiREJOKYWb5zzl9nu2YRCCFyzrF1XwmLCvexePM+Fm0OfN0dDInoKKN3h8Qjo4iBXZLp2ymRuGiFhIiELwVCPXHOsXnvoUBAVBtN7C0uByDGFwiJ7LTk4HWJ1vTumKCQEJGwoUBoQM45Cvcc+tpUU8Hmfew79P9Dom/HpEBABEcSvTskEhvdPO7yFZHwokBoZM45Nu0OhMSizXuPhERRSQUAsb4o+nYKTjcFp5x6d0gkxqeQEJGGpUAIA845NnxV/LWRxOLN+ygqDYZEdBT9OiUxsEsS2V1aM6BLMpkdEhQSIlKvFAhhqqrKsWF3MYsK9x65LrFky34OBEMiLhgSR65JpCXTKzWBaIWEiJwgBUIEqapyrPvq4NcuXC/ZvI+DZZUAxMdEkdUpiey01kdComdqAr4o87hyEYkECoQIV1nlWLfrIAWb91JQuJ+CzXtZsmU/xcGQaBHjo3/nr1+47lE9JGp574WINE8KhCaossqxdueBr93+unTLfg6VB0KiZWwgJM4vWs81d07GV1GOxcZCbq5CQaQZCzUQmv7SFU2IL8rI7JBIZodEJg5LA6Cisoo1Ow8GL1rvpWDzPvZ88CGUlWGuisrSMko/zKWlAkFE6qBAiHDRvij6dEykT8dELj0cEtkQNeZ1KsvKKIvyce2GVgydvoLrzupBcssYjysWkXAV0q0rZjbWzFaY2Wozu7WW4zeb2VIzW2RmuWaWXu3Yw2a2xMyWmdmTZmY1zp1qZotPvityWPSZZxD1SS6+++7lq3feJ+Xcs/n9jNWc+fAnPJm7iqKScq9LFJEwVOc1BDPzASuBMUAhMBe40jm3tFqbUcC/nXPFZnYjMNI5d7mZDQceAUYEm34G3Oacmxk87xLgUiDbOTegrmKb+zWEk7F0y34e/3glHy3dTuuWMdwwoieThqfTMlaDRJGmLtRrCKGMEE4FVjvn1jrnyoDXgIurN3DOzXDOFQc3vwTSDh8C4oFYIA6IAbYHC0wAbgbuC6EGOUlZnZN44Ro/U286g8FdW/PQv5Yz4uEZvDh7LSXBi9Ii0ryFEghdgE3VtguD+45mMjANwDk3B5gBbA0+pjvnlgXb3Qv8Fiiu7UkOM7PrzSzPzPJ27twZQrlyLNlprXn5u6fy1o059OmYyH3vL+PsR2bwlznrKa1QMIg0Z6EEQm3vfqp1nsnMrgL8BKaJMLNeQD8CI4YuwGgzG2Fmg4Fezrkpdb24c+5555zfOedPTU0NoVwJxbD0Nvzt2tN59brT6damJXe+u4TRj87itf9spLyyyuvyRMQDoQRCIdC12nYasKVmIzM7F7gdGOecKw3ungB86Zw74Jw7QGDkcDqQAwwzs/UEriv0NrOZJ9oJOXE5Pdvy+g05/Pl7p9IuMY5b3y7gnN/O4q38QiqrIuc9KiJy8kIJhLlAppl1N7NY4ApgavUGZjYEeI5AGOyodmgjcLaZRZtZDHA2sMw596xzrrNzLgM4E1jpnBt58t2RE2FmjOidyjs/GM5Lk/wkxkfzszcWMubxWUxduIUqBYNIs1BnIDjnKoCbgOnAMuB159wSM7vHzMYFmz0CJABvmNkCMzscGG8Ca4ACYCGw0Dn3Xn13QuqHmXFOvw7880dn8oerhhITFcX/vjqfC56czb8WbyOS3tUuIsdPS1fIUVVVOf5ZsJUnPl7J2p0HGdAliZvH9GZUn/bUeDuJiISx+rztVJqpqChj3KDOfPiTETx62SD2H6rgey/nMeGZL5i9aqdGDCJNjEYIErLyyirezC/kqdxVbNlXwqkZbbj5vN6c3qOt16WJyDFotVNpMKUVlfxj7iZ+/8lqdhSVckavttw8pg/D0lO8Lk1EaqFAkAZXUl7JX7/cwB9mrWHXgTJG9Unl5jF9GJiW7HVpIlKNAkEaTXFZBa98sYHnPl3D3uJyzsvqwE/H9KZfpySvSxMRFAjigaKScv742XpenL2WotIKLszuxE/PzaRX+0SvSxNp1hQI4pl9xeW8MHstf/p8HYfKKxk/uAv/e04mGe1aeV2aSLOkQBDP7T5YxnOz1vDKnPWUVzouHZrGj87pRVpKS69LE2lWFAgSNnYUlfDszDX87d8bcc5x+SlduWlUJh2T470uTaRZUCBI2Nm67xBPz1jNP+Zuwsz49mnduHFkT9onKhhEGpICQcLWpt3FPPXJKt6at5kYnzEpJ4Mbzu5Jm1axXpcm0iQpECTsrdt1kCdzV/HOgs20jPHxvTO7c+1ZPUhuEeN1aSJNigJBIsbqHUU8/vEq3l+0lcT4aK47qwffPSODxHgFg0h9UCBIxFm6ZT+Pf7ySj5Zup3XLGG4Y0ZNJw9NpGRvtdWkiEU2BIBFrUeFeHvtoJTNX7KRdQizfP7snV52eTnyMz+vSRCKSAkEiXv6GPTz20Qo+X/0VHZLiuGlUL755SlfiohUMIsdDgSBNxpdrv+KxD1fyn/W76dK6BT8a3YuJw9KI8enjPERCoUCQJsU5x+xVu/jtRytZuGkv6W1b8r+jMxk/pAu+KH16m8ix6BPTpEkxM0b0TuWdHwznpUl+EuKi+dkbCxnz+CzeW7iFqqrI+cNGJFxphCARyTnH9CXbePyjVazYXkTfjonc3W4/p21chI0aBTk5XpcoEjZCHSHofj6JSGbG2AGdOC+rI/8s2MqHL05h0M9vpqqqgqi4OCw3V6Egcpw0ZSQRLSrKGDeoM7/ruI+4qgp8VVW40jKYOdPr0kQijgJBmgTf6FFYXByVUT5Ko3zkZ2R7XZJIxNGUkTQNOTlYbi6VuZ9wx1dt+HB5NG9vLyKzgz6tTSRUGiFI05GTQ+wdt/Ozu75DXIyPya/ksedgmddViUQMBYI0OZ1bt+D5a4axbX8JN/4tn/LKKq9LEokICgRpkoZ2S+HBSwby5drd3DV1CZF0e7WIV3QNQZqsS4amsXL7Af4waw19OyZyTU6G1yWJhDWNEKRJ+/k3+nBuv/b8+r2lfLZql9fliIQ1BYI0ab4o44krhtArNYEf/C2fdbsOel2SSNhSIEiTlxAXzYuT/ET7opj8ylz2HSr3uiSRsBRSIJjZWDNbYWarzezWWo7fbGZLzWyRmeWaWXq1Yw+b2RIzW2ZmT1pASzN738yWB489WJ+dEqmpa5uWPPvtoWzaXcxNf59Hhe48EvkvdQaCmfmAp4HzgSzgSjPLqtFsPuB3zmUDbwIPB88dDpwBZAMDgFOAs4PnPOqc6wsMAc4ws/NPvjsiR3daj7bcN34As1ft4v4PlnldjkjYCWWEcCqw2jm31jlXBrwGXFy9gXNuhnOuOLj5JZB2+BAQD8QCcUAMsN05V+ycmxE8twyYV+0ckQZz+Snd+N4Z3fnT5+t59T8bvS5HJKyEEghdgE3VtguD+45mMjANwDk3B5gBbA0+pjvnvvanmZm1Bi4CckMvW+TE/fKCvozoncqd7yzmy7VfeV2OSNgIJRBq+ziqWt/lY2ZXAX7gkeB2L6Afgb/+uwCjzWxEtfbRwKvAk865tUd5zuvNLM/M8nbu3BlCuSLHFu2L4qkrh9CtbUtu/Gs+m3YX132SSDMQSiAUAl2rbacBW2o2MrNzgduBcc650uDuCcCXzrkDzrkDBEYOp1c77XlglXPuiaO9uHPueeec3znnT01NDaFckbolt4jhpUmnUOVg8itzKSrRnUcioQTCXCDTzLqbWSxwBTC1egMzGwI8RyAMdlQ7tBE428yizSyGwAXlZcFz7gOSgZ+cfDdEjl/3dq14+ltDWbPzID95bQGV+hhOaebqDATnXAVwEzCdwC/z151zS8zsHjMbF2z2CJAAvGFmC8zscGC8CawBCoCFwELn3HtmlkZgNJEFzAuec2299kwkBGdmtuOui7LIXb6Dh6cv97ocEU+FtJaRc+4D4IMa+35V7ftzj3JeJXBDLfsLqf3ahEijuyYng5Xbi3hu1lp6t09k4jDd8CbNk96pLALcdVF/cnq05ba3C8jfsMfrckQ8oUAQAWJ8UTzz7aF0ah3PDX/JY/PeQ16XJNLoFAgiQSmtYnlpkp/S8iqueyWP4rIKr0sSaVQKBJFqerVP5MlvDWH5tv387PWFVOnOI2lGFAgiNYzq055fXtCPaYu38UTuKq/LEWk0+sQ0kVpMPrM7K7YV8WTuKjLbJ3DRoM5elyTS4DRCEKmFmXHfhAGckpHCLW8sZFHhXq9LEmlwCgSRo4iL9vHsVcNolxDHdX/OY/v+Eq9LEmlQCgSRY2iXEMeLk/wUlVRw/Z/zKCmv9LokkQajQBCpQ79OSTxx+WAWbd7H/725COd055E0TQoEkRCc178jt5zXh6kLt/DMzDVelyPSIHSXkUiIfjCyJyu3F/HI9BX0TE1g7ICOXpckUq80QhAJkZnx0MRsBnVtzc2vL2Dplv1elyRSrxQIIschPsbHC1cPIyk+huv+nMfOotK6TxKJEAoEkePUPimeF67x89XBUr7/13xKK3TnkTQNCgSREzAwLZnfXjaY/A17uH3KYt15JE2CAkHkBF2Y3Ykfn5PJm/mFvDh7ndfliJw03WUkchJ+fE4mq3YU8Ztpy+jVPoFRfdt7XZLICdMIQeQkREUZj142iKxOSfzo1fms2l7kdUkiJ0yBIHKSWsZG88I1fuJjfEx+JY89B8u8LknkhCgQROpB59YteP6aYWzbX8KNf8unvLLK65JEjpsCQaSeDO2WwoOXDOTLtbu5a+oS3XkkEUcXlUXq0SVD01i5/QB/mLWGPh0SmTQ8w+uSREKmEYJIPfu/b/Th3H7tueefS5m9aqfX5YiETIEgUs+ioownrhhCr9QEfvi3eazdecDrkkRCokAQaQAJcdG8OMlPtC+Ka1/JY19xudclidRJgSDSQLq2acmz3x7Kpj3F3PTqPCp055GEOQWCSAM6rUdb7hs/gNmrdnH/B8u8LkfkmHSXkUgDu/yUbqzYdoA/fr6O3h0SufLUbl6XJFIrjRBEGsEvL+jLiN6p3PnOYr5c+5XX5YjUSoEg0giifVE8deUQurVtyY1/zWfjV8VelyTyXxQIIo0kuUUML006hSoH1/55LkUluvNIwktIgWBmY81shZmtNrNbazl+s5ktNbNFZpZrZunVjj1sZkvMbJmZPWlmFtw/zMwKgs95ZL9IU9a9XSue+fZQ1uw8yE9eW0BllZa3kPBRZyCYmQ94GjgfyAKuNLOsGs3mA37nXDbwJvBw8NzhwBlANjAAOAU4O3jOs8D1QGbwMfZkOyMSCc7o1Y67L8oid/kOHp6+3OtyRI4IZYRwKrDaObfWOVcGvAZcXL2Bc26Gc+7wpOiXQNrhQ0A8EAvEATHAdjPrBCQ55+a4wApgfwbGn3RvRCLE1TkZXHV6N56btZa38gu9LkcECC0QugCbqm0XBvcdzWRgGoBzbg4wA9gafEx3zi0Lnl/9f0FdzynS5Nx1UX9yerTltrcLyN+wx+tyREIKhNrm9mud+DSzqwA/8EhwuxfQj8CIoQsw2sxGHOdzXm9meWaWt3OnFgqTpiPGF8Uz3x5Kp9bx3PCXPDbvPeR1SdLMhRIIhUDXattpwJaajczsXOB2YJxzrjS4ewLwpXPugHPuAIGRw+nB50yrdnqtzwngnHveOed3zvlTU1NDKFckcqS0iuWlSX5Ky6u49pU8DpZWeF2SNGOhBMJcINPMuptZLHAFMLV6AzMbAjxHIAx2VDu0ETjbzKLNLIbABeVlzrmtQJGZnR68u+ga4N166I9IxOnVPpEnvzWEFdv287PXF1KlO4+Obc4ceOCBwFepV3UGgnOuArgJmA4sA153zi0xs3vMbFyw2SNAAvCGmS0ws8OB8SawBigAFgILnXPvBY/dCLwIrA62mVZPfRKJOKP6tOeXF/TjX0u28cTHK70uJ+xUVTlWbS/iw5emUD5yNO7OO+GccxQK9SyktYyccx8AH9TY96tq3597lPMqgRuOciyPwK2oIgJMPrM7K7cX8eQnq8nskMhFgzp7XZJnDpVVsrBwL/kb9hx57DtUzg/mvMPo8jLMVUFZGcycCTk5XpfbZGhxO5EwYWbcO34A63Yd5JY3FpLetiXZaa29LqtR7CgqIX/9HvI2BB5LNu+jIjh11jO1FWP7d2RYRgpnjIjH/ft1KsrKsJgYfCNHelt4E2OR9EHgfr/f5eXleV2GSIPadaCUi3//ORVVVUy96Uw6JMV7XVK9qqpyrNxRRN76wF/+eRt2s2l34A6ruOgoBqW1ZlhGCv70FIZ2SyGlVezXzi+f/Tl/efAV5nQbyGOPf5/E+BgvuhFRzCzfOeevs50CQST8LNu6n4nPfkFm+wT+cUMO8TE+r0s6YcVlFSzYuJe84NTPvI17KCoJ3E3VLiEOf3oKw9JTGJaRwoDOycRG132vy/yNe7jk2S+4+vR07rlYM891CTUQNGUkEob6dUriicsHc8Nf8/n5m4t48orBRMpyX9v2lZC3YfeREcDSrfuPrNnUu0MC/5PdGX96Cv6MFLq1aXlC/RrSLYXvDM/g5S/WM25QZ/wZbeq7G82SRggiYezpGat5ZPoKbjmvNzeNzvS6nP9SWeVYvm1/YOonGACH32AXHxPF4K6t8ae3YVhGCkO7ppDcsv6mdw6WVnDe458SHxPFBz8+i7joyB1FNTSNEESagB+M7Mmq7UU8+uFKerVPZOyAjp7Wc6C0gvkb//+dP/M37uVA8M10HZLi8Ke3YfKZ3RmWnkJW5yRifA23wn6ruGjunzCA7/xpLk9/spqbz+vTYK/VXCgQRMKYmfHgxGzWfVXMT/+xgK5tcujfObnRXn/z3kPkrd99ZASwfNt+qhyYQd+OSYwf0jkwAkhPIS2lRaNPa43s054JQ7rwzMw1XJDdib4dkxr19ZsaTRmJRIAd+0sY9/vPiTJ496YzSU2Mq/fXqKisYtnWosD8/4Y9zNuwh637SgBoGetjSLfWDEtvgz89hSHdWofN3T27D5Zx7mOz6NamJW/dOBxfVGRca2lMustIpIkpKNzHZc99Qf/Oyfz9utNOes58f0k58zfuJX99IAAWbNpLcVklAJ2T4xmW0ebIHUB9OyYS3YDTPyfr3QWb+fFrC/jV/2TxvTO7e11O2NE1BJEmZmBaMr+9bDA//Ps8bp+ymEcuzQ55isY5R+GeQ1+7+2fF9iKcgyiDrM5JfNPfNXD7Z3oKnVu3aODe1K9xgzrzzvzNPPrhCsZkdaBrm5ZelxSRFAgiEeTC7E6s3J7J73JX0btDAteP6Flru/LKKpZs2U/e+t3M2xiY/99RFFiEODEumiHpKZw/oBP+jBQGd21Nq7jI/lVgZtw3YSBjHpvF7e8s5pXvnhIxt+mGk8j+VyDSDP34nExW7SjigWnLGbp5Of71izhw+pnM7dTnyAhgYeFeSsqrAEhLacHwnm2PTAH17pDYJOfZu7Ruwf99ow93v7eUdxZsZsKQtLpPkq9RIIhEmKgo49HLBtEi7z/0v+anVFaW4/NF89QV97Ooaxb9OyfxrVPT8WcEpn+a2tIXx3J1TgZTF27hnveWMiIzlbYJ9X/xvSlTIIhEoJax0dyVuJPYynJ8roq4qkoe77CP9nd/gxaxzfcNWr6owG26Fz45m3v+uZTfXTHE65IiSvjeNiAix5R0/hh88XHg8xEVF0v6xAuadRgc1rtDIj8c1Yt3F2xhxvIddZ8gRygQRCJVTg7k5sK99wa+6nMBjrhxZE8y2ydw+5SCI++klropEEQiWU4O3HabwqCGuGgfD07MZuv+Eh6dvsLrciKGAkFEmqRh6SlMysnglTnryd+wx+tyIoICQUSarFu+0YdOSfH84q1FlFZUel1O2FMgiEiTlRAXzf2XDGT1jgM8M2ON1+WEPQWCiDRpo/q05+LBnXlm5mpWbi/yupywpkAQkSbvV/+TRUJcNL94a9GRT2+T/6ZAEJEmr21CHL+6KIv5G/fylznrvS4nbCkQRKRZGD+4C2f3TuXh6SuOfMynfJ0CQUSaBTPj/gkDALh9SgGR9FkwjUWBICLNRlpKS245rw8zV+xk6sItXpcTdhQIItKsTBqeweCurfn1e0vZfbDM63LCigJBRJoVX5Tx0MRsikrKufefS70uJ6woEESk2enTMZEbR/ZiyvzNzFyhFVEPUyCISLP0w1E96dU+gdunLOagVkQFFAgi0kzFRft4aOJAtuw7xKMfakVUCDEQzGysma0ws9Vmdmstx282s6VmtsjMcs0sPbh/lJktqPYoMbPxwWPnmNm84P7PzKxX/XZNROTYhqW34erT03n5i/XM26gVUesMBDPzAU8D5wNZwJVmllWj2XzA75zLBt4EHgZwzs1wzg12zg0GRgPFwIfBc54Fvh089nfgjnroj4jIcfn5N/rQMSmeW99aRFlFldfleCqUEcKpwGrn3FrnXBnwGnBx9QbBX/zFwc0vgbRanudSYFq1dg5ICn6fDOimYBFpdInxMdw3fgArtx/g2ZnNe0XUUAKhC7Cp2nZhcN/RTAam1bL/CuDVatvXAh+YWSFwNfBgCLWIiNS7c/p14KJBnfn9jFWsasYrooYSCFbLvlrf821mVwF+4JEa+zsBA4Hp1Xb/FLjAOZcG/Al47CjPeb2Z5ZlZ3s6dO0MoV0Tk+N11URat4qK59e0CqprpiqihBEIh0LXadhq1TO+Y2bnA7cA451xpjcPfBKY458qDbVOBQc65fweP/wMYXtuLO+eed875nXP+1NTUEMoVETl+7RLiuPPCLPI37OGv/97gdTmeCCUQ5gKZZtbdzGIJTP1Mrd7AzIYAzxEIg9re5XElX58u2gMkm1nv4PYYYNnxFi8iUp8uGdqFszLb8dC05Wxphiui1hkIzrkK4CYC0z3LgNedc0vM7B4zGxds9giQALwRvI30SGCYWQaBEcasGs95HfCWmS0kcA3h5/XSIxGRE2Rm/GbCQKoc3PHO4ma3IqpFUof9fr/Ly8vzugwRaeJenL2W+95fxpNXDmHcoM5el3PSzCzfOeevq53eqSwiUsN3z+jOoLRkfj11CXua0YqoCgQRkRp8UcaDE7PZd6ice99vPiuiKhBERGrRr1MSN47sydvzNvPpyuZxy7sCQUTkKH44qhc9UlvxyykFzWJFVAWCiMhRxMf4eGhiNoV7DvHYRyu9LqfBKRBERI7hlIw2XHV6N/70+ToWbNrrdTkNSoEgIlKHX4ztS/vEpr8iqgJBRKQOifEx3Dt+AMu3FfHcrKa7IqoCQUQkBGOyOnBhdiee+mQ1q3cc8LqcBqFAEBEJ0d0X9adFrI/b3l7UJFdEVSCIiIQoNTGOOy7sx9z1e/jbfzZ6XU69UyCIiByHS4elcWavwIqoW/c1rRVRFQgiIsfh8IqoFVVV3NnEVkRVIIiIHKdubVvyszF9+HjZDt4v2Op1OfVGgSAicgK+e0YG2WnJ3N2EVkRVIIiInIBoXxQPXpLN3uJy7v+gaXzgowJBROQEZXVO4oaze/BmfiGfrdrldTknTYEgInISfjQ6kx7tWnHblEUUl0X2iqgKBBGRkxAf4+OBSwayafchHo/wFVEVCCIiJ+m0Hm351mndeOmzdSyM4BVRFQgiIvXg1vP7kpoYxy/eWkR5ZWSuiKpAEBGpB0nxMdx7cWBF1Oc/Xet1OSdEgSAiUk/O69+RCwZ25He5q1izM/JWRFUgiIjUo7vH9Sc+Oorb3i6IuBVRFQgiIvWofWI8d1yYxX/W7ebVuZG1IqoCQUSknl3mT2N4z7Y8+MFytu0r8bqckCkQRETqmZnxwCUDKaus4s53I2dFVAWCiEgDSG/bipvH9OajpduZtnib1+WERIEgItJAJp/ZnQFdkvjVu0vYWxz+K6IqEEREGki0L4qHJmazp7iM30TAiqgKBBGRBtS/czLXj+jB63mFfL46vFdEVSCIiDSwH5+TSUbbltz2dgGHyiq9LueoQgoEMxtrZivMbLWZ3VrL8ZvNbKmZLTKzXDNLD+4fZWYLqj1KzGx88JiZ2f1mttLMlpnZ/9Zv10REwkNgRdRsNu4u5omPw3dF1DoDwcx8wNPA+UAWcKWZZdVoNh/wO+eygTeBhwGcczOcc4Odc4OB0UAx8GHwnO8AXYG+zrl+wGsn3x0RkfCU07MtV57alRdmr6WgcJ/X5dQqlBHCqcBq59xa51wZgV/cF1dvEPzFXxzc/BJIq+V5LgWmVWt3I3CPc64q+Bw7TqQDIiKR4tbz+9EuIXxXRA0lELoAm6ptFwb3Hc1kYFot+68AXq223RO43MzyzGyamWWGUIuISMRKbhHDPRcPYOnW/bwwO/xWRA0lEKyWfbW+7c7MrgL8wCM19ncCBgLTq+2OA0qcc37gBeCPR3nO64P8/DcLAAAGMElEQVShkbdz584QyhURCV9jB3RkbP+OPPHxKtbtOuh1OV8TSiAUEpjrPywN2FKzkZmdC9wOjHPOldY4/E1ginOuvMbzvhX8fgqQXduLO+eed875nXP+1NTUEMoVEQlv91zcn7joKG59a1FYrYgaSiDMBTLNrLuZxRKY+plavYGZDQGeIxAGtV0LuJKvTxcBvEPgQjPA2UD4XnoXEalH7ZPiuf2Cfvx73W7+kbep7hMaSZ2B4JyrAG4iMN2zDHjdObfEzO4xs3HBZo8ACcAbwdtLjwSGmWUQGGHMqvHUDwITzawAeAC49iT7IiISMS4/pSun92jDbz5Yxvb94bEiqkXKKnwAfr/f5eXleV2GiEi9WLfrIGOf+JSRfVJ57mp/g72OmeUHr9cek96pLCLike7tWvGTc3szfcl2/rV4q9flKBBERLx03Vnd6d85iTvfXcK+4vK6T2hACgQREQ8dXhF198EyHpjm7YqoCgQREY8N6JLMtWd157W5m/hijXcroioQRETCwE/O6U16cEXUknJvVkRVIIiIhIEWsT4emDCQDV8V88THqzypQYEgIhImhvdqx+X+wIqoizc3/oqoCgQRkTDyywv60aZVLL94axEVjbwiqgJBRCSMJLeM4Z5x/VmyZT8vfrauUV9bgSAiEmbGDujIeVkdePyjlaxvxBVRFQgiImHGzLh3/ABifVHc9nYBjbXEkAJBRCQMdUiK57YL+lE6+zMKbroV5sxp8NeMbvBXEBGRE3JFRSET/3EHvopy3J+ewnJzISenwV5PIwQRkTAV9eksYqsqiHZVUFYGM2c27Os16LOLiMiJGzkSi40Fny/wdeTIBn05TRmJiISrnBzIzQ2MDEaObNDpIlAgiIiEt5ycBg+CwzRlJCIigAJBRESCFAgiIgIoEEREJEiBICIigAJBRESCrLEWTaoPZrYT2HCCp7cDvPuwUm+oz82D+tz0nWx/051zqXU1iqhAOBlmluec83tdR2NSn5sH9bnpa6z+aspIREQABYKIiAQ1p0B43usCPKA+Nw/qc9PXKP1tNtcQRETk2JrTCEFERI6hyQSCmf3RzHaY2eJq+y4zsyVmVmVm/hrtbzOz1Wa2wsy+0fgVn5zj6a+ZtTWzGWZ2wMx+703FJ+84+zzGzPLNrCD4dbQ3VZ+c4+zzqWa2IPhYaGYTvKn65Bzv/+Xg8W7Bf9+3NG619eM4f84ZZnao2s/6D/VVR5MJBOBlYGyNfYuBS4BPq+80syzgCqB/8JxnzMzXCDXWp5cJsb9ACXAnEJH/Wap5mdD7vAu4yDk3EJgE/KXBq2sYLxN6nxcDfufc4OA5z5lZJC5x/zKh9/mwx4FpDVhTQ3uZ4+vzGufc4ODj+/VVRCT+Y6mVc+5TM8uosW8ZgJnVbH4x8JpzrhRYZ2argVOBhv8U63pyPP11zh0EPjOzXo1VX0M4zj7Pr7a5BIg3s7jgzzxiHGefi6ttxgMReYHwOP8vY2bjgbXAwUYor0Ecb58bSlMaIRyPLsCmatuFwX3SNE0E5kdaGJwIMzvNzJYABcD3nXMVXtfUkMysFfAL4Nde19LIupvZfDObZWZn1deTNpkRwnGqLXIj8q8pOTYz6w88BJzndS2NwTn3b6C/mfUDXjGzac65Eq/rakC/Bh53zh1ozL+kPbYV6Oac+8rMhgHvmFl/59z+k33i5hoIhUDXattpwBaPapEGYmZpwBTgGufcGq/raUzOuWVmdhAYAOR5XU8DOg241MweBloDVWZW4pyL2Jsn6hIc6ZYGv883szVAb+rh59xcA2Eq8HczewzoDGQC//G2JKlPZtYaeB+4zTn3udf1NAYz6w5scs5VmFk60AdY721VDcs5d2S6xMzuBg405TAAMLNUYLdzrtLMehD4/bW2Xp7cOdckHsCrBIZS5QRGAJOBCcHvS4HtwPRq7W8H1gArgPO9rr8R+rse2A0cCLbJ8roPDdln4A4CFxkXVHu097oPDdznqwlcQF8AzAPGe11/Q/e5xnl3A7d4XX8j/JwnBn/OC4M/54vqqw69U1lERIDme5eRiIjUoEAQERFAgSAiIkEKBBERARQIIiISpEAQERFAgSAiIkEKBBERAeD/Ad2ylezicFuoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "call_ms_16 = pd.read_csv('D:\\\\data\\\\options\\\\call_ms_16.csv',encoding='utf-8')\n",
    "\n",
    "price_ms_16=call_ms_16['最新价']\n",
    "k_ms_16=call_ms_16['执行']\n",
    "\n",
    "s0=112.35\n",
    "\n",
    "sigma_init=1\n",
    "sigma_16_newton=[]\n",
    "sigma_16_dichotomy=[]\n",
    "for i in range(call_ms_16.shape[0]):\n",
    "    sigma_16_newton.append(bsm_call_imp_vol_newton(s0,k_ms_16[i],t_16,rf,price_ms_16[i],sigma_init))\n",
    "    sigma_16_dichotomy.append(bsm_call_imp_vol_dichotomy(s0,k_ms_16[i],t_16,rf,price_ms_16[i]))\n",
    "print(sigma_16_newton)\n",
    "print(sigma_16_dichotomy)\n",
    "print('###############')\n",
    "plt.plot(k_ms_16,sigma_16_dichotomy,label='16-days',lw=1.5,)\n",
    "plt.plot(k_ms_16,sigma_16_dichotomy,'r.')\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.29000000639091406\n",
      "0.29\n",
      "0.15310000044057803\n",
      "0.1531\n",
      "0.002399991555676495\n",
      "0.0024\n"
     ]
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.36386969983929796, 0.36470650510148267, 0.3379696522663533, 0.308750828041821, 0.2901035211205926, 0.2770478308678258, 0.2678504434116783, 0.2795829812145216, 0.2815159982582912, 0.2892477867510631, 0.29487784902377207, 0.29826737962723826, 0.2737745536488281, 0.3638689603157307, 0.3647065048442553, 0.33796965226626546, 0.3087508280418194, 0.29010352112059345, 0.27704783086782625, 0.2678504434116786, 0.27958298121452063, 0.2815159982582914, 0.2892477867510613, 0.2948778490238011, 0.3001603015476018, 0.3079952199130588, 0.3638689603157307, 0.3647065048442553, 0.33796965226626546, 0.3087508280418194, 0.29010352112059345, 0.27704783086782625, 0.2678504434116786, 0.27958298121452063, 0.2815159982582914, 0.2892477867510613, 0.2948778490238011, 0.3001603015476018, 0.3079952199130588]\n",
      "###############\n",
      "[0.2663830055838007, 0.2623349770104382, 0.25886692974504427, 0.2551476532936361, 0.2538825758547821, 0.2511274813774724, 0.25311127192458355, 0.25282431711916775, 0.2533452730322892, 0.2663830055838, 0.2623349770104388, 0.25886692974504383, 0.25514765329363565, 0.25388257585478174, 0.25112748137747243, 0.25311127192458244, 0.25282431711916714, 0.25334527303228827, 0.2663830055838, 0.2623349770104388, 0.25886692974504383, 0.25514765329363565, 0.25388257585478174, 0.25112748137747243, 0.25311127192458244, 0.25282431711916714, 0.25334527303228827]\n",
      "###############\n",
      "[0.3052560705964372, 0.29001996299422883, 0.2942675861955737, 0.2691701404493931, 0.2708636644397628, 0.25993171220348904, 0.2610946078237797, 0.25271019230357544, 0.2523987108427236, 0.2475245983939655, 0.2483057133362069, 0.24353321749378337, 0.24308705976680645, 0.2385100817078438, 0.2406859629151303, 0.240161175751856, 0.30525607059637466, 0.2900199629942286, 0.29426758619557364, 0.26917014044939297, 0.2708636644397624, 0.25993171220348904, 0.2610946078237793, 0.2527101923035752, 0.2523987108427234, 0.24752459839396573, 0.2483057133362068, 0.24353321749378323, 0.24308705976680656, 0.23851008170784155, 0.2406859629151288, 0.24016117575185575, 0.30525607059637466, 0.2900199629942286, 0.29426758619557364, 0.26917014044939297, 0.2708636644397624, 0.25993171220348904, 0.2610946078237793, 0.2527101923035752, 0.2523987108427234, 0.24752459839396573, 0.2483057133362068, 0.24353321749378323, 0.24308705976680656, 0.23851008170784155, 0.2406859629151288, 0.24016117575185575]\n"
     ]
    }
   ],
   "source": [
    "for i in range(call_50etf_16.shape[0]):\n",
    "    sigma_16_newton.append(bsm_call_imp_vol(s0,k_16[i],t_16,rf,price_16[i],sigma_init))\n",
    "print(sigma_16_newton)\n",
    "print('###############')\n",
    "#\n",
    "\n",
    "for i in range(call_50etf_51.shape[0]):\n",
    "    sigma_51_newton.append(bsm_call_imp_vol(s0,k_51[i],t_51,rf,price_51[i],sigma_init))\n",
    "print(sigma_51_newton)\n",
    "print('###############')\n",
    "\n",
    "#\n",
    "\n",
    "for i in range(call_50etf_79.shape[0]):\n",
    "    sigma_79_newton.append(bsm_call_imp_vol(s0,k_79[i],t_79,rf,price_79[i],sigma_init))\n",
    "print(sigma_79_newton)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "x and y must have same first dimension, but have shapes (13,) and (26,)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-79-535a1826ce61>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mk_16\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0msigma_16_newton\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mlabel\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'16-days'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mlw\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1.5\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mk_16\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0msigma_16_newton\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'r.'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mk_51\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0msigma_51_newton\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mlabel\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'51-days'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mlw\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1.5\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mk_51\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0msigma_51_newton\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'g.'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mk_79\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0msigma_79_newton\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mlabel\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'79-days'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mlw\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1.5\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\matplotlib\\pyplot.py\u001b[0m in \u001b[0;36mplot\u001b[1;34m(*args, **kwargs)\u001b[0m\n\u001b[0;32m   3356\u001b[0m                       mplDeprecation)\n\u001b[0;32m   3357\u001b[0m     \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 3358\u001b[1;33m         \u001b[0mret\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0max\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   3359\u001b[0m     \u001b[1;32mfinally\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   3360\u001b[0m         \u001b[0max\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_hold\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mwashold\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\matplotlib\\__init__.py\u001b[0m in \u001b[0;36minner\u001b[1;34m(ax, *args, **kwargs)\u001b[0m\n\u001b[0;32m   1853\u001b[0m                         \u001b[1;34m\"the Matplotlib list!)\"\u001b[0m \u001b[1;33m%\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mlabel_namer\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfunc\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__name__\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1854\u001b[0m                         RuntimeWarning, stacklevel=2)\n\u001b[1;32m-> 1855\u001b[1;33m             \u001b[1;32mreturn\u001b[0m \u001b[0mfunc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0max\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1856\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1857\u001b[0m         inner.__doc__ = _add_data_doc(inner.__doc__,\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\matplotlib\\axes\\_axes.py\u001b[0m in \u001b[0;36mplot\u001b[1;34m(self, *args, **kwargs)\u001b[0m\n\u001b[0;32m   1525\u001b[0m         \u001b[0mkwargs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mcbook\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnormalize_kwargs\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0m_alias_map\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1526\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1527\u001b[1;33m         \u001b[1;32mfor\u001b[0m \u001b[0mline\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_get_lines\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1528\u001b[0m             \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0madd_line\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mline\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1529\u001b[0m             \u001b[0mlines\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mline\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\matplotlib\\axes\\_base.py\u001b[0m in \u001b[0;36m_grab_next_args\u001b[1;34m(self, *args, **kwargs)\u001b[0m\n\u001b[0;32m    404\u001b[0m                 \u001b[0mthis\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[0margs\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    405\u001b[0m                 \u001b[0margs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0margs\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 406\u001b[1;33m             \u001b[1;32mfor\u001b[0m \u001b[0mseg\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_plot_args\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mthis\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    407\u001b[0m                 \u001b[1;32myield\u001b[0m \u001b[0mseg\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    408\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\matplotlib\\axes\\_base.py\u001b[0m in \u001b[0;36m_plot_args\u001b[1;34m(self, tup, kwargs)\u001b[0m\n\u001b[0;32m    381\u001b[0m             \u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0my\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mindex_of\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtup\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    382\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 383\u001b[1;33m         \u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0my\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_xy_from_xy\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0my\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    384\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    385\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcommand\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;34m'plot'\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\matplotlib\\axes\\_base.py\u001b[0m in \u001b[0;36m_xy_from_xy\u001b[1;34m(self, x, y)\u001b[0m\n\u001b[0;32m    240\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m!=\u001b[0m \u001b[0my\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    241\u001b[0m             raise ValueError(\"x and y must have same first dimension, but \"\n\u001b[1;32m--> 242\u001b[1;33m                              \"have shapes {} and {}\".format(x.shape, y.shape))\n\u001b[0m\u001b[0;32m    243\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mndim\u001b[0m \u001b[1;33m>\u001b[0m \u001b[1;36m2\u001b[0m \u001b[1;32mor\u001b[0m \u001b[0my\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mndim\u001b[0m \u001b[1;33m>\u001b[0m \u001b[1;36m2\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    244\u001b[0m             raise ValueError(\"x and y can be no greater than 2-D, but have \"\n",
      "\u001b[1;31mValueError\u001b[0m: x and y must have same first dimension, but have shapes (13,) and (26,)"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAADYBJREFUeJzt3HGI33d9x/Hny8ROprWO5QRJou1YuhrKoO7oOoRZ0Y20fyT/FEmguEppwK0OZhE6HCr1rylDELJptolT0Fr9Qw+J5A9X6RAjudJZmpTALTpzROhZu/5TtGZ774/fT++4XHLf3v3uLt77+YDA7/v7fX6/e+fD3TO/fH/3+6WqkCRtf6/a6gEkSZvD4EtSEwZfkpow+JLUhMGXpCYMviQ1sWrwk3wuyXNJnrnC7Uny6SRzSZ5O8rbJjylJWq8hz/A/Dxy4yu13AfvGf44C/7T+sSRJk7Zq8KvqCeBnV1lyCPhCjZwC3pDkTZMaUJI0GTsn8Bi7gQtLjufH1/1k+cIkRxn9L4DXvva1f3TLLbdM4MtLUh9PPvnkT6tqai33nUTws8J1K35eQ1UdB44DTE9P1+zs7AS+vCT1keS/13rfSfyWzjywd8nxHuDiBB5XkjRBkwj+DPDe8W/r3AG8WFWXnc6RJG2tVU/pJPkycCewK8k88FHg1QBV9RngBHA3MAe8BLxvo4aVJK3dqsGvqiOr3F7AX01sIknShvCdtpLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDUxKPhJDiQ5l2QuycMr3P7mJI8neSrJ00nunvyokqT1WDX4SXYAx4C7gP3AkST7ly37O+CxqroNOAz846QHlSStz5Bn+LcDc1V1vqpeBh4FDi1bU8Drx5dvAC5ObkRJ0iQMCf5u4MKS4/nxdUt9DLg3yTxwAvjASg+U5GiS2SSzCwsLaxhXkrRWQ4KfFa6rZcdHgM9X1R7gbuCLSS577Ko6XlXTVTU9NTX1yqeVJK3ZkODPA3uXHO/h8lM29wOPAVTV94DXALsmMaAkaTKGBP80sC/JTUmuY/Si7MyyNT8G3gWQ5K2Mgu85G0m6hqwa/Kq6BDwInASeZfTbOGeSPJLk4HjZQ8ADSX4AfBm4r6qWn/aRJG2hnUMWVdUJRi/GLr3uI0sunwXePtnRJEmT5DttJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNDAp+kgNJziWZS/LwFda8J8nZJGeSfGmyY0qS1mvnaguS7ACOAX8GzAOnk8xU1dkla/YBfwu8vapeSPLGjRpYkrQ2Q57h3w7MVdX5qnoZeBQ4tGzNA8CxqnoBoKqem+yYkqT1GhL83cCFJcfz4+uWuhm4Ocl3k5xKcmClB0pyNMlsktmFhYW1TSxJWpMhwc8K19Wy453APuBO4AjwL0necNmdqo5X1XRVTU9NTb3SWSVJ6zAk+PPA3iXHe4CLK6z5RlX9sqp+CJxj9A+AJOkaMST4p4F9SW5Kch1wGJhZtubrwDsBkuxidIrn/CQHlSStz6rBr6pLwIPASeBZ4LGqOpPkkSQHx8tOAs8nOQs8Dnyoqp7fqKElSa9cqpafjt8c09PTNTs7uyVfW5J+UyV5sqqm13Jf32krSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSE4OCn+RAknNJ5pI8fJV19ySpJNOTG1GSNAmrBj/JDuAYcBewHziSZP8K664H/hr4/qSHlCSt35Bn+LcDc1V1vqpeBh4FDq2w7uPAJ4CfT3A+SdKEDAn+buDCkuP58XW/luQ2YG9VffNqD5TkaJLZJLMLCwuveFhJ0toNCX5WuK5+fWPyKuBTwEOrPVBVHa+q6aqanpqaGj6lJGndhgR/Hti75HgPcHHJ8fXArcB3kvwIuAOY8YVbSbq2DAn+aWBfkpuSXAccBmZ+dWNVvVhVu6rqxqq6ETgFHKyq2Q2ZWJK0JqsGv6ouAQ8CJ4Fngceq6kySR5Ic3OgBJUmTsXPIoqo6AZxYdt1HrrD2zvWPJUmaNN9pK0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqYlDwkxxIci7JXJKHV7j9g0nOJnk6ybeTvGXyo0qS1mPV4CfZARwD7gL2A0eS7F+27Clguqr+EPga8IlJDypJWp8hz/BvB+aq6nxVvQw8ChxauqCqHq+ql8aHp4A9kx1TkrReQ4K/G7iw5Hh+fN2V3A98a6UbkhxNMptkdmFhYfiUkqR1GxL8rHBdrbgwuReYBj650u1VdbyqpqtqempqaviUkqR12zlgzTywd8nxHuDi8kVJ3g18GHhHVf1iMuNJkiZlyDP808C+JDcluQ44DMwsXZDkNuCzwMGqem7yY0qS1mvV4FfVJeBB4CTwLPBYVZ1J8kiSg+NlnwReB3w1yX8mmbnCw0mStsiQUzpU1QngxLLrPrLk8rsnPJckacJ8p60kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNDAp+kgNJziWZS/LwCrf/VpKvjG//fpIbJz2oJGl9Vg1+kh3AMeAuYD9wJMn+ZcvuB16oqt8HPgX8/aQHlSStz5Bn+LcDc1V1vqpeBh4FDi1bcwj4t/HlrwHvSpLJjSlJWq+dA9bsBi4sOZ4H/vhKa6rqUpIXgd8Ffrp0UZKjwNHx4S+SPLOWobehXSzbq8bci0XuxSL3YtEfrPWOQ4K/0jP1WsMaquo4cBwgyWxVTQ/4+tuee7HIvVjkXixyLxYlmV3rfYec0pkH9i453gNcvNKaJDuBG4CfrXUoSdLkDQn+aWBfkpuSXAccBmaWrZkB/mJ8+R7g36vqsmf4kqSts+opnfE5+QeBk8AO4HNVdSbJI8BsVc0A/wp8Mckco2f2hwd87ePrmHu7cS8WuReL3ItF7sWiNe9FfCIuST34TltJasLgS1ITGx58P5Zh0YC9+GCSs0meTvLtJG/Zijk3w2p7sWTdPUkqybb9lbwhe5HkPePvjTNJvrTZM26WAT8jb07yeJKnxj8nd2/FnBstyeeSPHel9ypl5NPjfXo6ydsGPXBVbdgfRi/y/hfwe8B1wA+A/cvW/CXwmfHlw8BXNnKmrfozcC/eCfz2+PL7O+/FeN31wBPAKWB6q+fewu+LfcBTwO+Mj9+41XNv4V4cB94/vrwf+NFWz71Be/GnwNuAZ65w+93Atxi9B+oO4PtDHnejn+H7sQyLVt2Lqnq8ql4aH55i9J6H7WjI9wXAx4FPAD/fzOE22ZC9eAA4VlUvAFTVc5s842YZshcFvH58+QYuf0/QtlBVT3D19zIdAr5QI6eANyR502qPu9HBX+ljGXZfaU1VXQJ+9bEM282QvVjqfkb/gm9Hq+5FktuAvVX1zc0cbAsM+b64Gbg5yXeTnEpyYNOm21xD9uJjwL1J5oETwAc2Z7RrzivtCTDsoxXWY2Ify7ANDP57JrkXmAbesaETbZ2r7kWSVzH61NX7NmugLTTk+2Ino9M6dzL6X99/JLm1qv5ng2fbbEP24gjw+ar6hyR/wuj9P7dW1f9t/HjXlDV1c6Of4fuxDIuG7AVJ3g18GDhYVb/YpNk222p7cT1wK/CdJD9idI5yZpu+cDv0Z+QbVfXLqvohcI7RPwDbzZC9uB94DKCqvge8htEHq3UzqCfLbXTw/ViGRavuxfg0xmcZxX67nqeFVfaiql6sql1VdWNV3cjo9YyDVbXmD426hg35Gfk6oxf0SbKL0Sme85s65eYYshc/Bt4FkOStjIK/sKlTXhtmgPeOf1vnDuDFqvrJanfa0FM6tXEfy/AbZ+BefBJ4HfDV8evWP66qg1s29AYZuBctDNyLk8CfJzkL/C/woap6fuum3hgD9+Ih4J+T/A2jUxj3bccniEm+zOgU3q7x6xUfBV4NUFWfYfT6xd3AHPAS8L5Bj7sN90qStALfaStJTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ18f+GmWq6NWLIwgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(k_16,sigma_16_newton,label='16-days',lw=1.5,)\n",
    "plt.plot(k_16,sigma_16_newton,'r.')\n",
    "plt.plot(k_51,sigma_51_newton,label='51-days',lw=1.5)\n",
    "plt.plot(k_51,sigma_51_newton,'g.')\n",
    "plt.plot(k_79,sigma_79_newton,label='79-days',lw=1.5)\n",
    "plt.plot(k_79,sigma_79_newton,'k.')\n",
    "\n",
    "plt.grid(True)\n",
    "plt.xlabel('Strike price')\n",
    "plt.ylabel('Implied volatility')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
